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ABSTRACT 


The  location  times  of  thirty  measurement  pointe  for  a  space*' 
crafe  in  a  circular,  planar  earth  orbit  are  varied  to  minimize  a 
cost  function,  the  sum  of  the  squared  components  of  position  uncer** 
tainty,  at  a  pre-deterrained  target*  In  addition,  the  optimum  sched« 
ule  of  horizon  references  for  the  otar*elevation  measurement  to  be 
used  at  each  point  is  determined  with  respect  to  the  same  cost  func** 
tion*  A  steepest-dsecent  computer  program  was  written  to  perform 
the  optimization  in  each  case*  It  is  shown  that  the  measurement 
times  collect  into  four  clusters  from  a  nominal  schedule  in  which 
they  arc  equally  spaced*  A  cost  reduction  greater  than  807.  is  re» 
all zed*  The  horizon** selection  procedure  defines  certain  areas  along 
the  trajectory  where  one  or  the  other  horizon  is  preferred*  When 
carried  out  ainwltaneoualy  with  a  time  optimization,  this  procedure 
results  in  only  a  slight  improvement  over  the  case  where  a  single 
horizon  is  used  for  each  measurement* 
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CHAPTER  1 


INTRODUCTION 

The  optimization  of  a  celestial  measurement  schedule  for  a 
manned  space  mission  has  become  a  subject  of  great  importance  in 
the  last  few  years*  Walter  F*  Denham  and  Jason  L*  Speyer  of  Raytheon 
Company  considered  thla  problem  in  a  recent  report  (Reference  3)* 

They  oought  to  minimize  the  position  estimation  uncertainty  at  the 
terminal  point  of  a  fra©- fall  mission  by  comparing  various  sequences 
of  star  and  star-horizon  measurements*  A  steepest- deecent  numerical 
procedure  was  used  to  obtain  the  optimization*  The  authors'  results 
showed  a  10%  improvement  over  a  schedule  earlier  proposed  by  Richard 
H*  Battin,  of  the  MIT  Instrumentation  Lab*  In  this  and  similar 
studios,  the  locations  in  time  for  the  various  measurements  ware  held 
fixed  and  spaced  at  nearly  equal  intervals*  The  purpose  of  this  the¬ 
sis  is  to  investigate  the  behavior  of  such  a  nominal  time  schedule 
as  the  measurement  timaa  arc  changed  to  decrease  terminal  position 
uncertainty*  The  model  used  is  a  planar,  circular  earth  orbit  vherc 
the  target  point  appears  in  the  first  revolution*  Only  one  type  of 
measurement  Is  considered,  the  star-elevation  angle,  and  use  of  both 
horizons  is  investigated*  In  addition,  the  time  optimisation  prob¬ 
lem  is  coupled  vith  s  horizon-selection  procedure,  to  compare  with 
the  single-hjrlzon  rood 3* 
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It  it  expected  that  there  are  certain  preferred  measurement 
positions  along  the  circular  path*  If  the  times  of  the  measure* 
meats  Arc  free  to  chonge»  they  should  eventually  cluster  about  these 
points  in  order  to  effect  a  reduction  in  position  uncertainty*  A 
similar  hypothesis  con  be  stated  for  the  horizon  selection*  It  is 
probable  that  there  are  certain  areas  along  the  planar  trajectory 
where  It  would  be  more  beneficial  to  use  one  horizon  instead  of  the 
other*  The  stcepest^descent  procedure  will  be  used  to  determine 
the  optimum  schedules  in  both  cases * 


CHAPTER  2 


STATEMENT  OF  THE  PROBLEM 


Th«  objective  of  changing  th«  measurement  times  is  to  decrease 
tiis  position  uncertainty  ot  tht  target*  The  expression  for  the  co¬ 
variance  matrix  of  estimation  errors  is  developed  in  Appendix  B* 


-1 

E  -  A 
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(2-1) 


Where: 

Efi  •  estimation  error  covariance  matrix  at  target 

Eq  •  initial,  estimation  error  covariance  matrix 

C  -  state  transition  matrix  (Appendix  C) 

£  *  measurement  vector  (Appendix  A) 

2 

q-  •  variance  of  measurement  error 


Certain  assumptions  are  made  in  stating  the  problem  which  simplify 
the  elements  of  the  above  equation,  making  it  easier  to  manipulate* 
As  stated  ia  Chapter  1,  a  planar,  circular  earth  orbit  la 
assumed  for  the  spacecraft*  The  in-plane  navigstional  problem  can 
be  considered  alone  slnco,  as  is  shown  by  Stem  (Ref*  8)  end  others, 
the  in-end  out-of-pirn,©  error  propagations  are  uncoupled*  The  in¬ 
herent  simplicity  of  t.he  circular  o-Mt  ia  especially  obvious  in  the 
reduction  of  Stem*s  formula  for  the  transition  matrix  (Appendix  C) 
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Co  a  lass  complicated  form*  At  explained  In  Appendix  C,  Che  local 
vartlcal  co-ordinate  system  was  chosen  to  coincide  with  Stom9a 
aqua t Iona* 

It  vaa  decided  to  aelact  a  target  point  in  the  first  orbit  eo 
that  the  reaultant  time  changca  would  be  more  clearly  defined*  The 
entire  trajectory  is  Included  in  a  central  angle  of  290°t  and  the 
acro« angle  reference  la  arbitrary* 

The  stsr-elsv&tion  measurement  la  a  reasonable  selection  since 
it  hea  been  found  to  be  superior  in  the  vicinity  of  a  planet  (Chap*  % 
Ref*  1)*  Also-  ha  characteristic  vector  of  the  measurement,  devel¬ 
oped  in  Appendix  A,  t  ua  the  same  simple  form  at  all  points  in  the 
trajectory,  whan  expressed  in  local -vert leal  co-ordinates* 

Ae  implied  in  aquation  (2-1),  the  variance  of  the  measurement 
error  it  assumed  constant:  for  all  measurement;##  This  seems  reasonable 
enough,  since  the  type  of  measurement  la  the  same  each  timo  and  it  is 
always  taken  at  the  arms  altitude*  In  order  to  provide  ample  spac¬ 
ing  for  en  adequate  eirrole  of  measurement  points,  the  altltuda  chosen 
for  ths  problem  la  11 1 000  miles*  From  this  altitude,  an  optical  in¬ 
strument  con  be  expected  to  be  about  a  mile  in  error  in  discerning 
the  horizon*  Considering  an  error  of  about  *8  miles,  the  angular 
variation,  as  shown  in  Figure  2*1,  is  given  in  equation  (2-2)* 


*8 

ioV<n>*-<«2 


f)B  -  .554  tnr  <2»2) 

Expressed  In  arc'seoandc,  this  value  is  about  3.5  X  10"  aaconde. 
Hence,  the  variance  unod  in  the  problem,  assuming  the  mean  of  maseurt- 
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ment  errors  is  cero,  la  12  X  10~d  mc2* 


Figure  2*1 

Tho  qufntl  ties  for  the  Initial  estimation  orrore  are  chosen  to  be 
five  miles  and  ten  miles  per  second  In  position  and  velocity  in  each 
co-ordinate  directive  These  errors  are  aseumed  unco: ’related  so  that 
the  Initial  estimation  error  covariance  matrix  is  diagonal. 

Since  the  quantity  to  be  optimised  is  the  position  uncertainty. 
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only  the  first  two  diagonal  elements  of  the  4  X  4  Bg  matrix  ore  con¬ 
sidered*  A  convenient  way  to  write  this  cost  function  is  given  in 
equation  (2-3)* 


Cost 


tr 


Q  «* 


(2-3) 


A  more  sophisticated  cost  function  for  a  manned  mission  might  be  a 
weighted  average  of  the  target  position  and  velocity  errors*  such 
as  that  UBcd  by  Denham  and  Speyer*  This  would  imply  a  different  Q 
then  that  usod  ebova*  Another  possible  cost  function  is  the  deter¬ 
minant  of  tho  Ea  matrix*  describing  the  volume  of  tho  target  error 
ellipsoid* 

Stated  briefly*  tho  problem  is  to  find  the  time  schedule*  out 
of  ell  possible  schedules  of  thirty  roeasurernentc,  that  minimises  the 
cost  function  given  in  equation  (2-3)*  The  nominal  schedule  has 
thirty  racaeurament  points*  spaced  at  an  interval  of  about  900  sec¬ 
onds  in  tins*  between  central  angles  of  sere  and  290^*  Similarly* 
the  horlson-aclection  problem  seeks  to  find  she  sequence  of  horison 
references  which  minimises  cost*  There  is  a  choica  between  two 
references  at  each  point*  Tht  nominal  schedule  in  this  case  is  the 
use  of  tho  ♦’right”  horison*  opposite  to  the  direction  of  motion*  at 
each  point*  The  method  of  solution  In  each  case  is  the  steepest- 
descent  numerical  procedure*  which  is  the  subject  of  the  next  chap¬ 


ter* 


CHAPTER  3 


APPLICATION  OF  STEEPEST" DESCENT 


The  statpast- descent*  or  -ascent*  method  i  j  one  of  e  number 
of  numerical  techniques  developed  over  a  century  ago  by  Cauchy  and 
others  of  that  era*  The  advent  of  the  high-speed  computer  has 
brought  many  such  procedures  beck  to  life*  Largely  responsible  for 
the  revival  of  stsepcst-descent  are  Kelley  and  Bryson  who*  working 
independently)  recognised  its  superiority  in  certain  classes  of  prob¬ 
lems*  It  eliminated  much  of  the  guesswork  associated  with  other 
methods  by  assuming  a  non-optiraal,  nominal  solution)  and  proceeding 
to  the  optimum  by  a  ccrico  of  linear*  incremental  changes*  The  nom¬ 
inal  solution  need  only  bo  a  reasonable  first  guess  and  may  or  may 
not  satisfy  the  boundary  conditions* 

An  analogy*  credited  to  Brysoo*  illustrates  the  method  quite 
well*  A  hiker*  climbing  a  mountain  in  e  dense  fog*  will  climb  where 
the  slope  rises  the  sharpeot  to  minimi s«  his  time  of  ascent*  Because 
of  the  fog*  he  rcunt  relocate  the  direction  of  steepest  rise  at  regu¬ 
lar  intervale*  In  equation  form*  the  direction  in  which  he  climbs 
from  hi©  starting  point  1st 


r 

*r» 


i  4 
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where  s  is  tha  function  describing  the  hill*  The  horizontal  did* 
tance  moved  in  5  certain  direction  it  directly  proportional  to  the 
•lope  in  that  direction* 


Ax*  X2  • 


Ay  -  y2  -  yx  -  K*yi 


(3-2) 


The  linearising  assumption  is  that  the  total  vertical  dletance  climbed 
equals  the  sum  of  the  computed  vertical  distances  for  the  x  and  y 
directions* 

As-  s  A  x  +  Ay  (3*3) 

X1  7i 

Tha  climb  x  will  decide  before  he  starts  hot?  far  he  will  climb  ver» 
tically  before  ra-ac  cessing  tha  direction*  Hence#  A  a  is  a  known 
quantity* 

A  s  S  K  (t^)2  +  (Sy^)2  (3-4) 

The  constant  K#  which  governa  the  horizontal  di stance »  can  than  be 
determined* 


K  “  — r1 — 2  (3*5) 

<V  ♦  (*ri)2 

The  climber  p  edicts  that  his  new  oltituda#  when  he  has  arrived  at 
point  2#  will  be  ♦  A  z*  The  actual  altitude  will  normally  be  less 
than  this  straight  lino  extrapolation  of  slops*  After  determining 
tha  new  direction  of  steepest  ascent#  the  climber  repeats  the  pro* 
ctdure  until  finally#  the  actual  change  in  altitude  is  much  less 
than  he  predicted#  indicating  he  la  approaching  the  top  of  tha  moun* 
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tain* 

Two  disadvantages  of  the  ateepest-asccnt  technique  ere  brought 
out  in  the  analogy.  The  proper  step  sin,  As,  it  important  because 
the  climber  may  miss  a  better  path  if  he  climbs  too  far  in  any  direc¬ 
tion.  Unfortunately,  a  reasonable  step  sice  can  only  be  selected  by 
a  trial  and  error  process.  Also,  the  climber  may  venture  onto  an 
isolated  peak  and,  because  of  the  fog,  think  he  has  reached  the  top. 

A  fresh  start  with  ner*  initial  conditions  la  the  only  way  to  effec¬ 
tively  reduce  the  probability  of  converging  on  a  local  maximum. 

In  References  3  and  4,  Dry  son  has  outlined  the  mathematical 
approach  to  a  series  of  general  problems.  His  formulation  of  a  prob¬ 
lem  without  constraints  will  be  considered  here  since  it  is  somewhat 
similar  to  the  thesis  problem. 

A  nominal  spacecraft  trajectory  is  postulated,  which  is  de¬ 
scribed  by  the  following  eet  of  ordinary  differential  equations. 

dyi 

2^**  "  b,  t)  1  "  I#  2,  .  *  .,  n  (3-6) 


The  know  quantities  f^,  arc*  functions  of  the  independent  variable 
t,  the  dependent  variables  y4  and  the  driving,  or  control,  function, 
b(t).  The  cost,  s  function  of  the  depe  ^ent  variables,  is  increased 
(or  decreased)  by  varying  b(t).  Variations  about  this  nominal  tra¬ 
jectory  are  considered  end  it  it  assumed  that  they  can  be  accurately 


described  by  first-order  differentials  in  the  perturbation  equation. 


— -i.  i  y.  *  —v~  b  b 


(3-7) 


The  partial  derivatives  in  (3-7)  are  evaluated  along  the  nominal  trt 
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J factory •  Tho  dependent  variables  are  functions  of  time  so  that  (3*7) 
implied  a  set  of  n  linear  equations  with  variable  coefficients*  The 
l  terms  represent  a  small  variation  of  the  dependent  variables 
from  their  nominal  time  history*  A  set  of  equations  od joint  to  (3*7) 
is  defined  in  equation  (3*8)* 


(3-8) 


The  pertlal  derivative  in  (3*8)  is  the  negative  transpose  of  the  aim* 
liar  quantity  in  equation  (3*7)*  The  reason  for  the  adjoint  equation 
Is  made  clear  in  the  following  sequence! 

t  «■!& <Syi>*Syi^i)  -  £  Li^*b 

1-1  i«l 


*  £  £  <4  Q  S yj  -  4  «Ti>  (3.9) 

i-1  J-l  *yJ  *yl 

The  double  summation  term  in  (3*9)  equals  Eero  since  only  the  Indices 
differ*  The  left  hand  side  of  the  equation  is  equivalent  to  the  time 
derivative  of  5  y^* 


Z  LiSyi 


i-i 


(3-10) 


The  expression  which  relates  incremental  changes  In  the  control  vari¬ 
ables  to  the  resulting  changes  In  the  dependent  variables  is  obtained 
by  integrating  equation  (3-10)  over  the  flight  time* 
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n 


1  Li  s  y» 

1-1 


T  T 

r 


2  4  |Il  S 

i*i 


(3-11) 


L.  (t)  2  T  L.  IL. 
b  L  i  TT 


1-1 


2 -  Lt  $  yt 

l-l 


Lw(t)  S  b(t)dt 


Co  *o 


Th«  quantity  Lb(t),  defined  in  (3-11 )9  Is  the  influence  function  as¬ 
sociated  with  the  control  function9  b(t)*  The  definition  of  the  ad* 
joint  variable  Lj9  defined  in  equation  (3-8) 9  la  Juatiflad  by  this 
simple  expression*  Lj  la  a  known  function  of  the  nominal  trajectory 
and  lta  boundary  condition  is  a  function  of  the  cost9  which  is  usually 
determined  at  the  terminal  point  of  the  fl'ght* 


Lt(T) 


^  Coat 

b  yi 


t-T 


Cost 


Cost 


[y(t)J 


(3-12) 


The  objective  is  to  relate  changes  In  cost  to  changes  in  the  control 
function  by  the  use  of  the  adjoint  variable*  By  deflnitlon9  the  dif¬ 
ferential  cost  change  is  a  sum  of  partial  derivatives*  Using  equa¬ 
tion  <3-12)1 
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Cost 


t »T 


fa  a  *» 


’  n  1 

• 

£-  4  o  »» j 

i-1 

■I  t-T 

„  J 

(3-13) 


t-T 


Substituting  (3*13)  Into  (3*11) s 

T 

Scott  •  ^  Lb(t)  $  b(t)dt  ♦ 


z  h  sfl 

i«i 


(3*14) 


t-t. 


The  Adjoint  variables  Lj9  can  be  interprotad  as  ths  Influence  func* 
tlons  for  ths  initial  conditions  of  tha  dependant  variables*  Ths 
£  Cost  term  in  equation  014)  is  pre-selected*  For  a  given  value 
of  S  Coat*  it  is  desirable  to  require  the  smallest  possible  changes 
in  the  driving  function  and  initial  conditions  so  that  the  linear 
perturbation  equation  is  valid*  Stated  another  way*  the  problem  is 
to  minimise  the  effect  of  the  second-order  b  b(t)  and  £  y^  terms  for 
a  constant  cost  change*  The  summation  term  in  equation  014)  can 
be  rewritten  as  a  dot  product  to  simplify  the  mathematical 


"E  4  &  *i 

i-i 


ho  *  S*o 


(3-13) 


4 

*2 


where 
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If  the  1th  initicl  condition  It  specified,  f  is  scro  and  the  itl1 
tens  does  not  contribute  to  the  dot  product*  The  variational  calcu¬ 
lus  problem  can  then  be  stated  as  follows; 

T 

minimise!  J  •  j  |<)b(t)|2dt  ♦  c(  ^SXo  ^ 

(3-16) 

subject  toi  £  Cost  •  A  •  constant 


The  positive  constant  in  (3-16)  is  chosen  to  make  the  dimensions 
compatible  in  the  J  expression*  The  problem  can  be  rewritten  using 
Lagrange  multipliers; 

J*  *  J  ♦  TV  (SCost  •  A)  (3-17) 


where  J9  is  the  quantity  to  be  minimised*  Substituting  from  equa¬ 
tions  (3-14)  and  (3-15)i 

T 


J9 


1\  Lb(t)  $  b(t) 


♦ 


(3-18) 


to  *  S2o  -1U 

J9  can  be  divided  into  three  parts,  a  function  of  $  b(t),  a  function 
of  a  constant* 


■f  np  b(t)J  dt  ♦  F2C  £Xq) 


♦  constant 


(3-19) 


To  minimise  J9t 
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^Fi 

YsbftT 


IV  Lb(t)  +  2$  b(t)  -  0 


(3-20) 


*  ^  io  ♦  2*^  SXg  "  0 


0-21) 


The  second  derivatives  in  (3*20)  end  (3*21 )  ere  both  positive  snd  s 
minitsum  for  J*  is  assured*  These  equations  show  that  the  smallest 
changes  in  the  driving  function  and  the  initial  conditions  to  result 
in  a  given  cost  increment  are  changes  proportional  to  their  respec¬ 
tive  influence  coefficients* 


$  b(t)  ■ 

■  -  y  •  K 

(3-22) 

£2o  " 

*  tz1*  " 

(3-23) 

Thai  sign  of  the  constant  K  is  chosen  positive  or  negative  for  a  eta- 
sired  cost  increase  or  decrease*  Substituting  equations  (3-22)  and 
(3°23)  into  (3-14)  results  in  the  cost  expression  as  a  function  of  K* 


s 


Coat 


(3-24) 


Since  S  Cost  in  pre-selected  and  Lb(t)  end  are  known  functions 
daflned  by  equations  (3-11)  and  (3-15)  the  unknown  K  is  determined 
by  aquation  (3-25)* 


s 


Cost 


K 


(3-23) 
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The  influence  coefficients!  L0(t)  and  Lq*  determine  the  nature  of 
the  changes*  end  K  def  nninec  the  direction  and  magnitude* 

k^NEW  "  b^OLD  *  K  **b^  (3*26) 

^°NEW  "  &0LD  *  £  io  0*27) 

Ae  mentioned  previously*  only  the  unspecified  initial  conditions  ere 
available  for  change* 

The  general  procedure  can  be  summarised  as  follows: 

(I)  A  reasonable  first  estimate  of  b(t)  is  chosen,  according 
to  the  particular  problem* 

(II)  The  partial  derivatives  of  the  known  functions  fj  with 
respect  to  the  dependent  variables  and  the  control  func¬ 
tion  ere  avalueted  along  the  nominal  trajectory* 

(ill)  The  adjoint  variables  Lj  are  determined  from  equation 
(3-8)*  integrating  backward  over  the  nominal  trajectory 
with  equation  (3-12)  as  initial  conditioner  The  influ¬ 
ence  function  L^(t)  can  then  be  computed  from  equation 
(3-11). 

(iv)  An  arbitrary  cost  change  is  chosen,  depending  on  the 
nature  of  the  problem*  A  value  of  5  to  10%  might  be  e 
reasonable  initial  value  for  S  Coat  if  a  substantial 
overall  coat  change  is  anticipated*  K  is  then  deter¬ 
mined  from  equation  (3-23)* 

(v)  The  new  control  fv.  ictlon  is  found  from  equation  (3-26) 
and  the  new  initial  conditions  from  (3-27)*  Equations 
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(3-6)  are  then  integrated  to  obtain  the  new  trajectory 
and  the  process  ie  repeated* 

(vi)  The  ratio  of  the  predicted  cost  change  to  the  actual 

coat  change  will  Increase  as  the  optimum  is  approached* 
When  this  ratio  becomes  greater  than  about  S,  the  values 
of  &  Cost  and/or  K  should  be  decreased  to  reduce  step 
sise*  In  this  way,  the  optimum  can  be  approached  as 
closely  as  desired* 

The  essential  pert  of  this  general  formulation  is,  of  course9 
ralating  control  f vsict ion  changes  to  the  resultant  cost  change*  The 
adjoint  variables  were  necessary  to  obtain  such  an  expression  because 
a  direct  relation  between  £  Cost  and  f  b(t)  di  1  not  exist*  In  the 
thesis  problem,  the  estimation  error  covariance  matrices,  £k,  ere 
analogous  to  the  f^  in  the  general  formulation*  For  the  time  ael co¬ 
don  procedure,  the  times  of  the  measurement  points,  t^,  correspond 
to  the  driving  function  b(t>*  The  schedule  of  meaiuranent  vectors, 
la  the  driving  function  for  horizon  selection*  The  cost  func¬ 
tion,  explained  in  Chapter  2,  is  the  for  both  cases*  If  it  con 
be  expressed  as  an  explicit  function  of  tk  and,  for  the  other  case, 
the  adjoint  equations,  defined  in  (3-8),  will  not  be  needed*  In¬ 
stead,  the  influence  functions  for  both  cases  would  be  defined  by 
the  following  equations} 

Sco.t  -  Z  H  5*k 

k»l  k*»l 


(3-28) 
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Z  Cost 


Z  (^JiK  -  i  iSK 

k-1  k«l 


(3*29) 


Equations  (3*28)  and  (3*29)  are  both  analogous  to  aquations  (3-14) • 

The  Initial  condltionr  in  the  thesis  problacr.  art  specified  so  that 
the  term  corresponding  to  the  rightmost  term  in  (3-14)  is  xero* 

Also*  the  summations  ore  used  since  the  driving  functions*  unlike 
b(t)*  are  not  continuous  functions*  Using  the  cost  function  explained 
in  Chapter  2*  and  the  Efl  expression  developed  in  Appendix  B*  the  in¬ 
fluence  coefficients  for  the  two  cases  ere  derived  in  Appendices  0 
and  E* 


Time  selection: 


Lk  - 


~2  £  C*  E«QE«  -  ■  % 


* 


Horlson  selection: 


(D-13) 


h  *  -  £  C*  E.QE.(ci)’1  (E-17) 

<T2 

For  the  time  optimisation  problem*  the  N  state-transition 
matrices  are  evaluated  from  the  nominal  schedule  using  equation  (C-4)* 
Due  to  the  symplectlc  properties  of  C*  the  inverse  can  be  found  us¬ 
ing  the  elements  of  C«  Since  C  is  a  function  of  time*  the  determi¬ 
nation  of  its  derivative  is  straightforward*  These  results  are 
given  in  equations  (C-7)  and  (C«5)»  The  estimation  error  covariance 
fuatrix  is  then  computed  from  equation  (3-26)*  Using  these  quantities* 
and  the  measurement  vector  determined  in  Appendix  A*  the  time  selec¬ 
tion  influence  coefficient  is  obtained  from  oquatlon  (0-13)*  Since 
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a  decrease  in  coat  la  desired*  tha  time  change  at  each  point  la  op* 
poalta  in  sign  from  the  corresponding  influence  coefficient*  In* 
atead  of  specifying  a  specific  coat  change  as  in  Bryson's  formula* 
tioo*  it  ia  more  convenient  to  first  apecify  a  maximum  time  incre* 
ment«  If  tha  resultant  coat  change  la  too  email*  the  time  increment 
can  be  increased  until  a  change  greater  then  JZ  le  obtained*  of  course* 
ae  the  optimum  le  approached*  a  smaller  percentage  change  ia  required* 
The  time  changes  ere  scaled  according  to  the  else  of  their  respective 
Influence  coefficients  ee  in  equation  (3*22)*  so  that  the  measurement 
point  having  the  greatest  effect  on  coat  ia  changed  the  most*  The 
change  follow  the  direction  of  steepeat  descent  in  an  N*diraeasional 
•pace*  The  intricacies  of  this  procedure  are  clarified  by  the  flow 
charts  In  the  next  chapter* 

The  horison- select  ion  procedure  ia  similar  up  to  a  point* 

Since  it  is  carried  out  simultaneously  with  the  time-optimisation 
procedure*  Efl  and  Cak  oust  be  re* evaluated  after  each  iteration*  There 
ia  not  an  much  control  in  this  problem  however*  since  there  ere  only 
two  possible  valued  for  £  at  each  point*  The  N  individual  elements 
of  equation  (3*29/  must  be  examined  to  determine  the  incremental  coot 
changes*  If  e  proposed  horizon  change  rosultn  in  a  decrease  in  coat* 

the  change  ia  made*  If  not*  the  original  horison  la  retained*  With 

0 

so  little  control*  it  le  possible  that  the  proposed  change  violates 
tho  assumption  of  llnaarity  utilised  in  the  perturbation  equation* 

This  problem  is  discussed  in  Chapter  5*  Tho  procedure  for  horison 
selection  is  also  illustrated  in  Chapter  4* 


CHAPTER  4 


COMPUTER  SOLUTION 


The  specifics  of  the  computer  program  used  to  implement  the 
theory  developed  in  Chapters  2  and  3  and  Append lcea  A  thru  E  ere 
covered  in  Appendix  F#  However*  it  will  be  useful  to  understand  how 
the  problem  solution  is  carried  out#  The  flow  charts  in  Figures  4#l» 
3  will  help  in  understanding  the  methods  used# 

Figure  4#1  gives  tho  flow  chart  for  Block  One  (no  horlson 
change)*  Here  the  now  measurement  time  schedule  is  oomputed  using 
the  same  measurement  vector*  either  left  or  right  horlson* 

The  input  data  needed  is  oovered  in  Appendix  F#  From  equa* 
tion  (A*14)»  the  measurement  vectors  can  be  computed* 


&(RIC!IT) 


&(LEFT) 


(4*1) 


(4*2) 
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Using  these  vectors  end  the  stata  transition  matrix  equation 
(C<>4)  the  target  estimation  error  covariance  matrix,  equation  (B-26) 
can  be  computed* 

As  shovn  in  Appendix  B,  the  symplectio  property  of  the  state 
transition  matrix  aUovs  one  to  simply  compute  the  matrix  inverse 
and  matrix  transpose  inverse  by  rearranging  the  elements*  Equation 
(C*»7)  is  computed  in  a  simple  subroutine* 

The  matrix  Q  in  covered  in  Appendix  D*  Nov,  the  cost  computed 
will  be  the  followings 

Cost  •  tr  |qE^J  (D*»2) 

The  cost  will  be  designated  the  old  cost,  oc,  when  the  computation 

A 

uses  a  measurement  time  schedule  which  is  either  the  Initial  one  or 
a  result  of  a  previous  iteration* 

Using  the  present  target  estimation  error  covariance  matrix 
and  equation  (D-13)  the  influence  coefficients  are  computed*  The 
logic  used  will  change  tha  measurement  time  schedule  by  an  amount 
depending  on  the  influence  coefficient  having  the  largest  magnitude* 
Now,  if  we  define  a  scale  factor,  sf,  ass 

mt  m  I  maximum  influence  coefficient! 
maximum  time  Increment 

or 

■f  -  l£Cml  (4»3) 

A  *m 

then  the  new  measurement  time  depends  on  the  old  measurement  time 
end  the  velue  of  tha  influence  coefficient  at  the  old  measurement 


tilkJ#  Novi 
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new  time  »  old  time  •  Influence  coefficient 

scale  fee cor 


or 


(4-4) 


cad  this  procedure  is  applied  to  all  the  meaaurenvmt  times*  Obviously 
the  time  with  tha  influence  coefficient  having  the  largest  magnitude 
will  have  he  greatest  change#  Also  the  sign  of  the  influence  coef- 
ficient  will  determine  which  way  the  measurement  time  will  move* 

After  all  the  measurement  times  have  been  changed,  a  new  tar* 
got  estimation  error  covariance  matrix  and  a  new  coet*  ncf  can  ba 
computed*  Tha  actual  cost  change*  acc,  is  clearly s 


acc  -  oc  «  nc  (4*5) 
and  thie  number  should  be  positive*  A  predicted  coat  change*  pcc, 
can  be  defined  aoj 


30 

pcc  -  2T  EC,  (At,) 
1-1  1 


30 

EC,  (t,o  -  t,n) 


(4-6) 


*5nd  then  it  is  compared  to  acc*  The  percentage  coat  change  la  them 


Pc  -  a*  (4-7) 

When  pc  . s  positive,  it  la  than  compared  to  soma  minimum  de¬ 
sired  percent  change*  rape  If  pc  is  less  than  mpe*  then  tha  maximum 
time  increment  is  multiplied  by  the  ratio  of  mpe  to  pc* 

When  dc  is  saro*  the  maximum  time  increment  it  cut  in  half* 
wiisn  pc  ie  negative*  the  new  maximum  time  increment  is  changtd  by 
an  amount  depending  on  hov  much  negative  it  ia*  The  logic  le  then 
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if  pc  +  rape  is  zero,  then: 

-  Aj”  (4-tO 

if  pc  +  rape  is  greater  than  zero  then: 

"  A*»<1  "  ^5  >  «-» 

end  if  pc  4-  rape  is  less  than  zero  then: 

A  *•  -  A  «»■  *  (4-10) 

Having  changed  the  maximum  time  increment,  it  may  be  used  to 
repeat  the  measurement  time  change  procedure  shown  previously*  The 
same  influence  coefficients  are  used  since  thoy  were  calculated  be* 
fore  any  measurement  time  change  was  made*  The  same  procedure  of 
selecting  a  scale  factor,  and  of  then  changing  cadi  measurement 
time  is  carried  out*  This  procedure  may  continue  until  pc  la  equal 
to  or  greater  then  mpe*  Any  additional  changes  made  will  be  added 
to  those  previously  made*  Since  tho  logic  does  not  return  to  the 
measurement  time  schedule  used  to  compute  the  influence  coefficients, 
the  optimum  may  be  mlesed?  much  as  the  hill  climber  in  Chapter  3 
related  a  better  path  by  climbing  too  far  in  ono  direction*  The  ini* 
tiel  time  increment  may  hove  forced  the  optimum  over  the  top  of  the 
hill  and  any  further  change  in  this  time  lncroment  will  merely  cause 
chmgca  that  will  put  the  optimum  further  ova?  the  top*  Zhs  program 
may  then  run  into  trouble*  When  this  happens,  the  best  procedure 
to  follow  is  to  return  to  the  measurement  time  schedule  computed  be* 
fore  the  trouble  was  encountered  and  reduce  cither  the  maximum  time 
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increment  or  tha  minimum  desired  percent  change* 

Therefore,  the  program  control  it  either  on  the  percent 
change  or  the  maximum  tine  increment*  Another  possible  program 
control  could  be  the  ratio  of  pcc  to  acc* 

The  above  procedures  are  shown  better  in  Figures  4*2  and 
4*3*  Subroutine  CHECK,  Figure  4*2,  changes  only  one  time*  This 
routine  is  used  right  after  the  influence  coefficients  are  computed* 
The  result  shows  that  the  coat  does  decrease  by  changing  the  measure* 
rcent  time,  and  that  the  Influence  coefficients  ere  correct* 

Subroutine  LOGIC,  Figure  4*3,  computes  the  entire  new  measure* 
rant  time  schedule*  Also  this  routine  limit*  the  new  timet  to  the 
end  conditional 

0  -  tj  -  fa  (4-11) 

n 

v^ere  Fa  is  tht  final  angle  of  the  orbit* 

This  whole  procedure  can  be  repeated  any  number  of  time** 

But  act  mentioned  before,  after  several  Iterations,  it  may  be  1  ropes* 
slble  to  achieve  e  given  rope*  Increasing  the  maximum  time  increment 
may  piece  the  optimum  ovor  the  hill  top*  Then  e  new  rape  must  be 
chosen  and  this  can  only  be  don*  in  a  heuristic  manner*  However, 
the  results  after  a  s »t  of  iterations  halp  determine  what  sics  steps 
must  be  made  to  bring  the  measurement  tiro*  schedule  to  an  optimum* 
Figure  4*4  shove  the  flow  chert  for  the  second  block*  Her* 
a  nr.#  measurement  tinn  schedule  and  a  new  measurement  horlson  schedule 
are  computed*  The  mein  difference  lies  in  the  feet  that  at  each  mens* 
urement  time  tha  measurement  vector  is  different,  and  the  results  of 
the  program  ere  op  tier  jm  measurement  time  end  horiaon  echetkilce* 
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Except  for  using  a  different  measurement  vector  at  each  ti  ae® 
the  measurement  time  schedule  optimization  ie  the  acne  as  in  Block 
One«  Also  since  linearity  ie  assumed  superposition  holds,  and  the 
measurement  time  and  horizon  schedule  optimisation  can  be  carried  out 
independently* 

The  measurement  time  schedule  optimisation  is  completed  first 
and  then  using  equation  (E-l 7 >  the  horizon  influence  vectors  ere  com* 
puted. 

Using  equations  (E«20)  the  change  in  coat  it  computed*  Only 
if  tills  change  in  cost  is  negative,  will  the  total  cost  be  reduced 
by  changing  tha  measurement  vector* 

To  avoid  going  outside  tha  linear  range  by  making  the  change 
in  cost  too  large,  only  one  measurement  vector  will  be  changed  at  a 
time*  This  will  keep  the  change  in  cost  small*  Therefore,  only 
the  measurement  time  having  the  negative  change  in  cost  with  ths 
largest  magnitude  will  have  its  measurement  vector  changed* 

After  this  particular  measurement  vector  is  changed,  the  new 
targot  estimation  error  covariance  matrix  end  a  now  coat  are  computed* 
Aa  beforei 

acc  »  oc  -  nc  (4-3) 


and  the  predicted  change  in  cost  ia  defined  by  equation  (4**12)t 


where] 

pcc  •  -  5  Cost 

(4-12) 

S  Cost  -  - 

<r  2 

^  ^ok  'a^a  ^ak  ^ 

(E-l  6) 

The  ratio  of  pcc  to  acc  shows  how  the  change  effected  the  cost* 
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The  new  coat  can  be  compared  to  the  old  cost  to  see  if  change  did 
decrease  the  cost* 

Subroutine  JUMP*  Figure  4*5  ihovs  the  logic  used  in  changing 
the  measurement  vector* 
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a 

acc 

ACk 

c 

a 

ECj 

ECM 

Fa 

Sficl 

$  &lr 
IH^ 


mpc 

n 

nc 

NI 

oc 

pc 

pcc 

r 

s£ 


TABLE  4-1 

SYMBOLS  FOR  FLOVJ  CHARTS 

ratio  of  actual  cost  change  to  old  cost 
actual  cost  change 

cost  change  due  to  change  in  measurement  vector 
estimation  error  covariance  matrix  at  target 
influence  coefficient  for  measurement  time  change 
largest  value  for  ECi 
final  angle  of  the  orbit 

measurement  vector  change  from  r.ght  to  left  horizon 
meacurement  vector  change  from  left  to  right  horizon 
horizon  flag 

-  0  using  left  horizon 

•  1  uaing  right  horizon 

maximum  desired  percent  coat  change 
mean  angular  motion 
new  cost 

number  of  iterations 
old  coat 

percentage  cost  change 
predicted  cost  change 

ratio  of  predicted  cost  change  to  actual  cost  change 

scale  factor 

new  measurement  time 


old  measurement  time 
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maximum  time  increment 
actual  time  increment 
influence  vector  • 

Subecriptai 

i9  k  measurement  times 

I  one  particular  measurement  tiroo 

N  mtaauremen:  time  having  the  a  -  with  the  largest 

magnitude 


& 

A 

— “k 


CHAPTER  5 


RESULTS  AND  CONCLUSIONS 

As  mentioned  in  Chapter  2,  the  hypothesis  associated  with 
the  time  optimisation  problem  is  thit*  in  a  schedule  of  thirty 
equally  spaced  measurement  points*  there  are  a  certain  number  of 
preferred  positions*  Mers jremants  made  at  or  near  these  positions 
should  result  in  a  lover  cost*  the  sum  of  the  squares  of  terminal 
position  uncertainty*  than  measurements  made  at  other  points  along 
the  trajectory*  If  the  times  of  the  various  points  are  allotted  to 
change  to  effect  a  cost  decrease*  they  should  cluster  about  the  pre« 
f erred*  or  optimum*  points*  The  method  of  steepest  descant  is  par* 
tlcularly  applicable  to  this  type  of  problem  since  the  relative  else 
of  the  thirty  Influence  coefficients  indicates  the  sensitivity  of 
their  corresponding  points*  The  time  changes  are  proportional  and 
of  opposite  sign  from  their  respective  coefficients*  so  that  a  rela¬ 
tively  large  value  for  Lk  indicates  that  a  substantial  tin#  change 
should  be  mane  in  a  specific  direction*  If  the  optimum  points  are 
vellodeflncd*  thalr  position  on  the  trajectory  should  not  change  ap¬ 
preciably  os  the  times  of  ths  measurements  are  changed*  Therefore* 
a  plot  of  the  influence  coefficients  as  s  function  of  the  correspond¬ 
ing  central  angles  for  each  case  should  serve  to  locate  them*  These 
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influence  function*  ere  plotted  in  Figure*  5*1  *nd  5*2  for  the  cost 
values  listed  In  Table  5**1* 


NOMINA! 

CHANGE  -  25% 

CHANGE  -  501 

RIGHT 

.6405  x  109 

*4686  x  109 

.3221  x  109 

LEFT 

1. 1628  x  109 

*8432  x  109 

.5732  x  109 

COST  VALUES  -  FT2 

TABLE  5-1 

Tie  "RIGHT'  horlson  is  that  opposite  to  the  direction  of  motion*  as 
defined  in  Appendix  A*  It  is  evident  from  Table  3**1  that*  for  a 
single  horizon  reference*  the  right  horizon  is  preferable*  The  an* 
gular  ('Isperalon  of  the  measurement  points  corresponding  to  the  above 
coat  talues  are  more  clearly  shown  in  the  polar  plote*  Figures  5*3 
through  5*7*  The  general  configurations  of  the  Influence  functions 
In  figures  5*1  and  5*2  remain  the  same*  even  after  substantial  changes 
in  cost*  The  Increase  in  amplitude  indicates  that  the  times  are 
driven  harder  toward  the  optimum  an  the  optimum  la  approached*  The 
arrows  in  both  figurea  Indicate  the  direction  of  time  change*  and 
serve  to  define  the  circled  stable  points*  It  appears  that  four  dus- 
tor/;  should  result*  tvo  at  the  end  points  and  two  In  the  middle*  The 
cluster  locations  predicted  from  Figures  5*1  and  5*2  are  given  In 


Table  5*2* 
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i»t 

2nd 

3rd 

4th 

RIGHT 

0° 

58.5° 

170° 

290° 

LEFT 

0° 

70° 

176° 

290? 

CLUSTER  PREDICTIONS 
TABLE  5-2 


For  a  total  coat  change  greater  than  50%,  the  dispersion  of  the 
points  is  not  sufficient  to  accurately  identify  the  sero- crossings 
in  an  influence  function  plot*  Subsequent  Iterations  were  carried 
out,  periodically  decreasing  the  maximum  time  Increment  and  required 
percent  change,  until  the  clusters  were  clearly  defined*  The  angular 
dispersion  for  a  cost  change  of  about  75%  is  shown  in  Figures  5-8 
and  5-9*  At  this  point,  it  is  obvious  that  there  will  be  only  four 
clusters*  The  Influence  coefficients  at  this  stage  tend  to  rive  a 
number  of  the  times  beyond  the  end  points*  As  noted  in  Chapter  4, 
the  measurement  positions  arc  constrained  once  they  reach  0°  end 
290~'»  and  the  lar^e  end-point  influence  coefficients  are  ignored  in 
computing  the  predicted  coet  change  for  each  iteration* 

Accurate  identification  of  the  position  of  the  dusters  was 
not  possible  until  after  several  iterations  requiring  a  *1%  cost 
decresso  or  loss*  The  size  of  the  clusters  cannot  ba  predicted  since, 
in  the  early  iterations,  the  program  drove  tha  time  locations  quits 
hard  until  a  substantial  cost  decrease  was  realised*  There  were 
osveral  Instances  of  points  "jumping”  from  one  duster  to  another* 

It  is  reasonable  to  assume  that  a  tighter  tolerance  on  the  maximum 
time  increment  would  result  in  different  cluster  si  sea*  The  loose 
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to l «r taco  tma  used  to  shortaa  the  convergence  time. 

Whoa  the  program  had  changed  the  time#  os  much  os  possible* 
the  most  likely  cluster  positions  vote  chosen*  and  all  points  were 
assigned  one  of  theso  four  time  values*  The  cost  function  values 
resulting  from  the  selected  solution  show  only  a  slight  decrease  from 
the  oooputer  solution*  The  results  are  listed  in  Table  5-3* 


COST 

COMPUTER 

COST 

SELECTED 

X  CHANGE 

RIGHT 

*1108  *  109 

*1108  x  109 

82.77. 

LEFT 

*2162  x  1G9 

*2161  x  109 

81.47* 

FINAL  COST 

VALUES 

TABLE 

5-3 

The  final  angular  positions  ond  numbers  of  included  points  for  the 
clusters  ore  given  in  Table  5-4*  along  with  the  positions  predicted 
from  Figures  5-1  end  5-2* 

RIGHT  LEFT 


Included 

Points 

Angle 

Pred. 

Included 

Points 

Angle 

Pred 

1st 

O 

a# 

0° 

0° 

3 

0° 

0° 

2nd 

11 

69.3° 

58.5° 

11 

75.4° 

70° 

3rd 

7 

208.1° 

170° 

10 

197.2° 

176° 

4th 

10 

290° 

290° 

6 

290° 

290° 

FINAL  CLUSTER  POSITIONS 
TABLE  5-4 


The  angular  position  of  the  clustors  is  more  clearly  shown  in  Figures 
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5-10  end  5-11*  The  influence  fimction  plot*  provide  *  fair  predict 
tion  o*  the  number  end  position  of  the  clusters  even  before  any  cost 
reduction  is  obtained* 

The  selected  cluster  positions  can  be  Justified  only  if  the 
influence  coefficients  for  these  points  approach  zero,  indicating 
that  there  are  no  further  change*  to  be  made*  The  influence  coeffi¬ 
cient*  for  th*  end  point*  are  *till  quite  large  but,  a*  noted,  they 
tend  to  drive  the  time*  beyond  the  con*traint*.  The  sign  of  the  co¬ 
efficients  1*  positive  at  0°  and  negative  at  290°,  so  that  the  times 
are  being  driven  in  th*  proper  direction*  The  value*  of  the  coeffi¬ 
cient*  for  the  middle  two  cluster*  are  oompa red  with  the  values  for 
point*  close  to  th*  cluster  position  in  Table  5-5* 


ON 

CLUSTER 

OFF 

CLUSTER 

ANGULAR 

DIFFERENCE 

RIGHT 

2nd 

♦  4*40 

-  68.0 

2.5° 

3rd 

-  6.36 

♦  69.3 

3.4° 

LEFT 

2nd 

-21.32 

-173.5 

3.9° 

3rd 

♦  5.59 

•  22.0 

4.6° 

INFLUENCE  COEFFICIENTS  -  FT2/ SEC 
TABLE  5-3 

It  ia  evident  that  the  optimum  positions  have  been  closely  approxi¬ 
mated* 

Ihe  cost  reductions  in  each  case  will  be  more  meaningful  if 
compared  in  terras  of  position  uncertainty  in  the  radial  and  tangen¬ 
tial  directions.  As  noted  in  Chapter  2,  the  initial  estimation  error 
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was  chosen  to  b«  about  five  miles  in  each  direction*  The  correspond* 
log  values  after  making  thirty  measurements  for  the  nominal  and  opti* 
mum  solutions  are  compared  in  Table  5*6* 


NOMINAL 

OPTIMUM 

%  IMPROVEMENT 

RIGHT 

RADIAL 

2*53 

1*33 

47.5% 

TANGENTIAL 

4*08 

1*49 

63.5% 

lept 

RADIAL 

2*98 

1*50 

49*6% 

TANGENTIAL 

5*74 

2*34 

59*2% 

POSITION  UNCERTAINTY  •  MILES 
TABLE  5-6 

Changing  the  times  of  the  measurement  points  results  in  a  significant 
cost  reduction  in  both  cnaos  and  the  right  horizon  reference  gives 
the  best  results* 

The  purpose  of  the  her 1 son-selection  procedure  was  to  Inves¬ 
tigate  the  posaibllity  or  e  cost  reduction  by  providing  e  choice  of 
two  references  st  each  measurement  point*  Using  the  measurement  vec¬ 
tors  defined  in  Appendix  A  for  the  right  and  left  horizons,  the  £  & 
vectors  were  defined  in  equations  (£-18)  end  (E-19)*  The  time  opti¬ 
mization  problem  described  in  the  first  part  of  this  chapter  was  car¬ 
ried  out  first*  Since  it  wee  evident  after  the  first  run  that  the 
right  horizon  reference  would  result  in  a  lower  cost  value,  a  reason¬ 
able  rnilnal  schedule  for  the  horizon  selection  was  to  use  this  hori¬ 
zon  «t  each  point*  The  idee  was  to  switch  to  the  left  horizon  where 
the  steepest  descent  procedure  predicted  a  decrease*  Since  the 
switching  problem  was  paired  with  a  timing  schedule  optimization,  it 
was  anticipated  that  the  horizon  schedule  would  not  stabilize  until 


the  optimum  time  schedule  was  approached*  Stated  another  way,  the 
portioi  t  of  the  trajectory  which  preferred  one  horizon  over  another 
were  expected  to  be  a  function  of  the  central  angle  only. 

Aa  explained  in  Appendix  S,  the  horizon  aelection  differed 
fro*;  the  time  optimization  in  that  there  waa  no  control  ovar  the  etep 
size*  The  measurement  vector  could  not  be  driven  in  a  direction  to 
effect  a  coat  decrease  since  the  two  values  cf  were  predetermined* 
If  the  cost  change  predicted  from  a  proposed  horizon  change  was  nega* 
tlvc,  the  switch  was  made*  If  not,  the  original  horizon  waa  retained* 
The  problem  does  not  have  the  continuous  nature  of  the  time  optimize* 
tion  and  the  lack  of  step  size  control  caused  trouble*  Early  results 
using  the  scheme  descrloed  above  did  not  provide  accurate  predictions 
of  cost  change.  When  the  program  changed  the  horizon  at  all  points 
whore  a  cost  decrease  was  predicted,  the  resultant  coat  valua  waa 
greater  than  before*  The  Influence  vectors  were  correct,  so  it 
seemed  beat  to  change  only  one  horizon  at  a  time  before  re-evaluating 
the  vectors*  The  problem  peraiated,  howavtr,  and  at  that  point  the 
atap  size  waa  investigated*  It  was  found  that  the  right  and  left 
horizon  vectors  were  separated  by  an  angle  of  149°  at  11,000  miles 
aa  shown  in  Figure  5*12* 
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COMPARISON  OF  S &  VECTOR  AT  DIFFERENT  ALTITUDES 


11,000 

nlle 

orbit 


If  the  altitude  was  reduced  to  its  lowest  practical  limit,  100 
miles,  :he  angle  is  reduced  to  12°;  but,  since  the  vectors  art 
inversely  proportional  to  the  altitude,  it  was  felt  that  the  step 
size  would  still  be  too  large*  The  alternative  solution  was  to  re¬ 
define  the  "left"  horizon  vector,  using  the  negative  of  the  £  (LEFT) 
defined  in  Appendix  E*  The  physical  meaning  of  this  change  is  that 
the  star-elevation  angle  nould  be  treasured  in  a  counter-clockwise, 
rather  than  clockwise,  direction*  An  examination  of  Figure  1  in 
Appendix  A  shows  that  this  is  true*  The  time  optimization  procedure 
is  not  affected  by  this  change  since,  in  the  expressions  for  EA, 

and  S  Cost  (Equations  (D-10),  (D-13),  (D«»l)),  the  measurement  vcc- 

T 

tors  appear  in  the  form  ^  •  Therefore  only  the  square  of  the  ole- 
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tr.cn ts  is  critical  and  changing  tho  sign  of  the  measurement  vector 
docs  not  weaken  the  comparison  with  the  straight  tlroe-optirai  cation 
problem*  Ac  shown  in  Figure  5-13*  the  step  sice  was  coiaiderably  de¬ 
creased*  The  step  sice  could  be  reduced  more*  if  accessary*  by  in¬ 
creasing  the  altitude* 


As  a  further  precaution*  only  one  horizon  change  was  made  before 
re-evalunting  the  Influence  vectors*  The  vectors  replacing 
those  in  Appendix  E  are  given  in  equations  (5*1)  and  (5-2)* 
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S\  (LEFT  TO  RIGHT) 


S  %  (RIGHT  TO  LEFT) 


2r, 


*<*2  -  r  *)* 


2rt 


2^ 


s(z2  -  r£2) 


0 

0 


(5*1) 


(5-2) 


The  procedure  described  in  Chapters  3  end  4  was  carried  out* 
with  more  encouraging  results*  The  predicted  incremental  cost  changes 
are  plotted  in  Figure  5*14  as  a  function  of  the  central  angle*  Since 
the  nominal  schedule  uses  the  right  horizon  at  all  points,  the  areas 
of  negative  cost  predict  a  favorable  change  to  the  left  horizon*  It 
is  evident  that  the  points  preferring  one  horizon  or  the  other  are 
not  scattered  randomly  over  the  trajectory  but  lie  together  in  cer¬ 
tain  well  defined  areas* 

Since  the  time  and  horizon  selection  procedures  are  indepen* 
dent,  the  shape  of  Figure  5-14  should  not  be  affected  by  changing 
measurement  times*  The  angular  limits  corresponding  to  either  hori* 
zon  reference  can  be  predicted  from  the  figure*  These  predictions 
are  listed  in  Table  5*7* 
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RIGHT  HORIZON  LEFT  HORIZON 

16°  to  109°  d°  to  16° 

I 

212°  to  290°  109°  to  212° 

AREA  PREDICTIONS 
TABLE  5-7 

A  similar  line  of  reasoning  applies  to  the  time-optimization 
clusters*  Changing  tho  horizon  from  right  to  left  should  not  affect 
the  number  of  clusters  although  their  positions  may  be  ollghtly  al¬ 
tered*  The  predictions  for  the  cluster  positions  are  given  in  Table 
5”8  along  vlth  the  predicted  horizon  obtained  from  Table  5-7* 


ANCLE 

HORIZON 

lot 

0° 

LEFT 

2nd 

72° 

RIGHT 

3rd 

202° 

LEFT 

4th 

290° 

RICHT 

CLUSTER  PREDICTIONS 
TABLE  5*8 

As  expected,  the  effect  of  two  opt *mi cation  procedures  is 
to  provide  for  more  rapid  convergence*  The  horizon  changes  made 
along  vlth  tho  respective  central  angles  are  listed  in  order  of 
their  occurrence  in  Table  5*9*  For  each  iteration,  the  proposed 
horizon  change  which  results  in  the  greatest  cost  decrease  is  the 
only  change  made*  The  numbers  associated  vlth  the  polnt6  are  the 
identification  numbers  in  the  program*  After  several  iterations, 
these  numbers  lose  their  meaning  since  the  points  any  pass  each 
other  on  the  vay  to  the  optimum* 
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POINT  NO* 

CHANGE 

ANGLE  AT  CHANGE 

19 

RIGHT  TO  LEFT 

163.5° 

18 

« 

137.0° 

17 

w 

148.5° 

20 

tt 

174. 0^ 

16 

140. (P 

21 

H 

182.4° 

15 

ti 

135.6° 

1 

M 

8.2° 

14 

n 

122.2° 

22 

190.8° 

13 

n 

112.1° 

13 

LEFT  TO  RIGHT 

112.0° 

2 

RIGHT  TO  LEFT 

17.2° 

22 

LEFT  TO  RIGHT 

240.9° 

21 

LEFT  TO  RIGHT 

HORIZON  CHANGES 

TABLE  5-9 

232.9° 

The  predictlone  in  Tabic  5-7  were  quite  accurate  for  several 
iteration r  At  one  point  however,  while  seeking  an  overall  cost 
reduction  of  5Z  or  greater t  the  program  made  changes  which  were  ob¬ 
viously  outside  the  linear  range*  A  number  of  the  angles  were 
changed  by  30°  or  more*  A  cost  reduction  was  realized  from  these 
new  values  but  the  large  changes,  in  effect,  altered  the  nature  of 
the  problem*  If  such  violations  of  linearity  were  not  allowed  the 
new  values  would  be  arrived  at  from  a  different  nominal  schedule* 

This  "new"  nominal  schedule  would  probably  result  in  different  zero- 
crossings  in  Figure  5-14  and  hc.co  different  predictions  In  Table 
5-7*  This  line  of  reasoning  seeks  to  explain  the  apparent  discrepan¬ 
cies  in  the  last  two  entries  of  Table  5-9* 

As  the  clusters  become  more  clearly  refined.,  the  horizon 
selection  stabilizes  since  there  is  no  further  movement  across  the 
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boundaries*  The  result  of  the  time  optimisation  is  given  In  Table 
5-10* 


INCLUDED  ANGLE  HORIZON 

POINTS 


1st 

2 

0° 

LEFT 

2nd 

H 

67° 

RIGHT 

3rd 

6 

209° 

LEFT 

4th 

11 

290° 

RIGHT 

FINAL  CLUSTERS 
TABLE  5*10 

The  table  shows  that  the  horizon  selection  procedure  does  not  change 
appreciably  the  strength  and  position  of  the  clusters*  Comparing 
Tables  5^8  snd  5*10  shows  the  accurscy  of  the  predictions*  The  die* 
persion  of  the  points*  at  different  stages  in  the  optimisation*  is 
shown  in  the  polar  plots*  Figures  5*15  through  5*18*  The  nominal 
positions  are  the  seme  as  show  in  Figure  5*3*  Note  that  the  Im¬ 
provement  in  overall  percent  change  is  %  over  the  BRIGHT*’  case 
in  Table  5*3*  which  has  the  tone  initial  conditions*  The  sizes 
of  the  time ’selection  Influence  coefficients  for  points  at  and  near 
the  cluster  positions  are  compared  in  Table  5*11*  in  order  to  Justify 
the  final  position  of  the  dusters* 


ON 

OFF 

ANGULAR 

CLUSTER 

CLUSTER 

DIFFERENCE 

1st 

15.0 

33.0 

o 

o 

• 

2nd 

92.8 

872.9 

.02° 

INFLUENCE  COEFFICIENTS  -  F?2/SEC 


TABLE  5-11 
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The  radial  and  tangential  components  of  the  final  position  uncer* 
talnty  are  listed  in  Table  5-12*  In  these  more  familiar  units,  the 
improvement  over  the  previous  time  optimisation  is  more  obvious  * 


NOMINAL 

OPTIMUM 

7.  IMPROVEMENT 

RADIAL 

2.53 

1.05 

50.2* 

TANGENTIAL 

lt.08 

1.56 

65.316 

POSITION 

UNCERTAINTY  -  MILES 

TABLE  5*12 


* 


rr 


'  ran  inn  naHHaran 


rr  ft* 


time  change  takes  on  a  different  meaning,  dependant  on  where  It  oc¬ 
curs  along  the  trajectory* 

The  time-optimization  problem  can  bo  coupled  to  a  more  sophis¬ 
ticated  measurement- selection  procedure  as  an  extension  for  the 
second  portion  of  the  thesis*  As  discussed  in  Chapter  3,  the  steep- 
est-descent  technique  is  not  always  a  rel labia  method  in  this  dis¬ 
continuous  type  of  problem*  The  step  size  from  one  possible  measure¬ 
ment  to  another  must  be  small  enough  so  that  the  linearizing  assump¬ 
tions  are  valid  or  false  predictions  will  result*  Denham  and  Speyer 
were  conscious  of  these  limitations  and  eteepcst-dcscent  work ad  well 
in  their  study* 

Different  cost  functions  can  be  used  to  determi  .e  their  effect 
on  the  clusters*  Some  possible  schemes  are  mentioned  in  Chapter  2, 
such  as  a  weighted  average  of  position  and  velocity  uncertainty* 
Instead  of  making  many  measurements  in  a  shovt  space  of  time,  as  is 
implied  by  the  cluster?  of  points,  it  might  be  more  reasonable  to 
track  the  angular  elevation  of  a  star  for  a  certain  amount  of  time* 
This  would  ellroi'  ate  froouent  changes  of  spacecraft  altitude  and  re¬ 
sult  in  a  significant  fuel  oavlng*  A  cost  function  could  be  con¬ 
trived  which  would  compare  the  effect  of  tracking  on  terminal  posi¬ 
tion  uncertainty  to  the  effect  of  making  discretes  measurements* 

It  is  interesting  to  speculate  on  the  physical  reasons  behind 
the  clucters*  It  might  be  argued  that  a  series  of  measurements  early 
in  the  flight  would  tend  to  reduce  the  effec-  of  the  initial  estima¬ 
tion  error  by  Improving  the  estimate  before  it  propagates  too  far* 
Likewise,  the  measurements  near  the  terminal  point  would  tend  to 


67 


reduce  the  position  uncertainty  of  that  point*  Ths  position  of  ths 
middle  two  clusters  is  somewhat  of  a  mystery,  however*  A  loo,  the 
strength  of  the  clusters,  measured  by  the  number  of  included  points, 
and  their  convergence  times  are  by  no  means  identical*  Tho  values 
may  be  random  in  nature,  but  it  is  more  likely  that  there  are  physi¬ 
cal  causes* 

The  results  obtained  are  a  strong  argument  for  tho  practicality 
of  changing  the  measurement  times,  in  spite  of  the  simplicity  of  the 
problem*  It  is  reasonable  to  expect  that  a  more  complicated  model 
would  produce  similar  results* 


APPENDIX  A 


MEASUREMENT  VECTOR 


The  measurement  vectorc  for  a  variety  of  techniques  are  da* 
vcloped  in  Reference  1*  For  this  problem,  tho  atar«elavation  measure* 
ment  was  chosen  because  of  its  superior  accuracy  in  o  planetary  orbit* 

The  local  vertical  co-ordinate  system  is  used  to  coincide  with  the 
system  adopted  by  Stern  in  developing  his  transition  matrix* 

As  shown  in  Oattin,  each  type  of  measurement  is  characterised 
by  o  measurement  vector*  The  deviation  in  the  quontity  to  be  measured 
relates  to  the  position  deviation  by  the  following  formula,  vhsre  h 
is  this  vector* 

$  q  «  h  •  J  r  (A*l) 

The  procedure  for  determining  h  is  the  oama  for  all  measure¬ 
ments;  the  equation  defining  the  quantity  to  bo  measured  is  developed 
and,  from  it,  the  perturbation  equation*  In  this  case,  the  measured 
quantity  is  the  angle  from  the  planet  hori2on  to  tho  line  of  sight 
to  a  star,  oc  shown  in  Figure  1*  Equation  (A-2)  defines  the  measured 
angle  M* 

n  •  z  «  2  coe(M*B)  (A-2) 

where  n  la  e  unit  vector  in  the  direction  of  the  star*  Writing  the 
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V3 


STAR 


Figure  1 


STAR  ELEVATION  ANGLE 


perturbation  equation) 


n  •  g  z  «  -2  sin(M+B)(  £m+  £  B)  ♦  £  a  CO«(M+B)  (A-3) 

A  vector  expression  for  the  scaler  82  vill  be  useful* 

z  •  z  o 

S  z  •  ^  •  £  ac  •  2s  £z 

z  •  £  £ 

2  •  8 2  (A-4) 

Substituting  (A-4)  into  (A»3)j 

z  •  £  a 

= - -  cos(M«-B)  -  2  oin(M+B)(SM+  $B)  °  n  •  g  z  (A«5) 

2 

From  Figure  29  it  is  obvious  that  &Z"  -  Jr* 


r 

Planet 
Center 

1  Position 

Actual  position 


FIGURE  2 


r  +  £r  •  r* 

z  •  £r  «  «  j*  ♦  £  £  by  definition 

Sr  «  -  £  2 

Rewriting  (A-5)  in  terms  of  Sr: 

z  ain(M+B)(  6m+  £b)  ■  £  n  -  m  cob(M+B)  j  ♦  £r  (A-6) 

where  a  is  the  unit  vector  in  the  direction  of  planet  center*  From 
Figure  1: 

cos(M+B)  «  m  •  n  the  projection  of  n  on  m 

Therefore,  the  bracketed  quantity  in  (A-6)  relates  tue  two 
legs  of  a  right  triangle*  Figure  3  shows  that  the  magnitude  of  the 
resultant  lag  is  sln(M+B>* 


1  -  cos2(M+B)  «  sin2(M+B) 


FIGURE  3 


Rewriting  (A-6)  in  term?  of  the  vector  at 


£  M  +  &  B  «  £  *>  l 

t 


An  expression  for  $3  can  be  derived  from  Figure  !• 


•in  B  « 


D 

2  r 


Writing  the  perturbation  of  (A-8): 


cos  B  S  B  ** 


**  £  Z 

2z2 


Substituting  (A-4)  and  the  results  of  Figure  2  in  (A-9)i 


S  B 


■  D  £  •  S  £ 


2z^cos  B 


S  B 


Cm 


23^008  B 


Sr 


Substituting  (A”10)  into  (A<*7): 


g  M  o  £  •  S  r  ^  Dm  *  g  r 
2  2z2coa  B 


From  Figure  Is 


oin  B  «  ? 


2z 


J  M  ■  *,(n  •»  tan  Bn)  •  $  r 
z  ~  ~  “ 


(A -7) 


(A-8) 


CA-9) 


(A*10) 


(A-ll) 


It  is  shown  in  Figure  3  that  a  is  perpendicular  to  m*  Figure  4  shows 
that  the.  bracketed  quantity  in  equation  (A~ll)  defines  the  third  leg 
of  a  right  triangle* 


FIGURE  4 

tan  Q  •  •  tan  B 

a 

Q  •  B 

b  la  perpendicular  to  d  9  the  unit 
vector  In  tho  direction  of  the  planet 
edge* 

Determining  the  magnitude  of  b  i 

1  ♦  tan^B  ■  b*  •  sec^B 
Defining  the  unit  vector  £  3 

b 

£m  — S.,. 

sec  B 

Equation  (A— 11 )  can  be  rewritten! 

Sh  -  *• 


z  cos  B 


(A.12) 


Referring  to  aquation  (A-l),  the  vector  which  characterises 
the  star  elevation  measurement  is: 


h 


8  COS  B 


(A-13) 


As  noted  in  Figure  4,  £  is  perpendicular  to  the  planet  edge 
end  is  therefore  independent  of  the  measured  angle  M*  Hence*  the 
expression  for  the  deviation  in  M  does  not  contain  M  explicitly* 

The  h  vector  can  be  written  directly  from  the  orbital  geometry* 

The  assumption  of  a  circular  orbit  at  a  known  altitude  serves 
to  completely  determine  the  measurement  vector*  The  expression  is 
derived  in  Figure  5*  in  the  local  vertical  co-ordinate  system,  rsz* 


r 


FIGURE  5 


LOCAL  VERTICAL  CO-ORDINATE  SYSTEM 


£  "  -  cos  Bs  +  sin  Br 

sin  B  *  rE 
z 

cos  8  -  («2-rE2)* 

z 

From  (A-13). 

h  .  •  cos  Bji  -f  sin  Br 
Z  COS  B 


h 


“  I  8  +  i  tsn  Br 
z  “  Z  — 


h 


rE 

z(»  -rE2)^ 


r 


(A-14) 


If  the  left  horizon  is  used,  the  s«  component  changer*  sign* 

In  order  that  it  be  compatible  with  the  4X4  state  transition 
matrix,  the  measurement  vector  is  modified  to  the  4* dimensional  vec« 
tor  £»  vhoee  first  two  elements  are  identical  with  h* 


£ 


Equation  (A»l)  then  becomes: 


(A-15) 


&.  •  s* 


(A-16) 


where  £x  is  the  complete  state  deviation  vector* 


APPENDIX  B 


ERROR  CORRELATION  MATRIX 


The  error  vector  e  ma  be  defined  os  the  difference  between 
Che  estimated  end  true  velues  of  the  state  deviation  vector  at  any 
point  in  the  trajectory* 

e  -  -  5x  (B»  I) 

-n  v  -n  °  -ti 

The  covariance  matrix  of  the  error  vector  io  of  interest  sinco  it 
describes  the  uncertainty  in  the  estimated  position  and  velocity 
deviation  of  the  spacecraft* 

En  ”  SnSnT  <B*2> 

Battin  (Reference  2)  develops  a  recursion  formula  for  the  error  cor* 
relation  matrix  which  will  be  presented  in  this  section*  It  is  then 
shown  that  the  recuroion  formula  is  equivalent  to  a  more  convenient 
representation  using  the  Inverse  of  the  covariance  matrix*  It  is 
this  latter  form,  referred  to  target  co-ordinates,  which  la  utilized 
in  the  computer  program*  Definition  of  the  symbols  to  be  used  is 
necessary  before  proceeding*  The  following  superscripts,  applying 
to  the  quantities  Sx,  b  q  and  c,  ore  defined! 

obsarvation 
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9  extrapolation 
^  estimation 

The  lack  of  a  superscript  on  these  symbols  indicates  their  trie 
values*  The  extrapolated  vector  is  that  resulting  from  the  same 
vector  at  a  previous  time,  if  no  nev  measurements  intervene* 


&  %  “  cn,n-l  S  5n-l 


(B-3) 


Equation  (B-S)  defines  the  state- transition  matrix,  C*  The  extrapo¬ 
lated  error  matrix  is  encountered  in  the  subsequent  derivation  and 
is  defined  as  follows: 

•J  -  $  k  ■  S  24  “  Cn,n-1  §  i»-l  *  Cn,n-1  &  2o~l 

e*  «  C  .  e  . 

— nfn»’l  — n-1 


e*e,T 

-n-n 


Cn,n-i  ^n-1 


1 


(8-4) 


An  equation  similar  to  (B-4)  transfers  any  matrix  from  ito  local  co¬ 
ordinate  system  to  a  reference  system,  target  co-ordinates)  in  this 

case* 


C  EC 
an  n  an 


(B-5) 


The  seme  relation  must  hold  true  for  the  extrapolated  matrices* 


Substituting  from  equations  (8-4)  and  (B-5): 

E*®  -  C  C  n  .  E  .  C  I  , 

rt  an  n,n-l  n-1  n,n-l  an 


(3-6) 
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^fl,n-i  ^-1  ^a,n-i 


(B-7) 


Thus,  the  estimation  error  covariance  matrix  at  point  n,  just  prior 
to  the  nth  measurement,  equals  the  matrix  just  after  the  (n-i)**1 
measurement  when  they  are  referred  to  the  same  co-ordinates*  This 
is  a  clearer  definition  of  the  extrapolated  covariance  matrix  than 
that  which  equation  (B  4)  provides* 

A  linear  estimn:e  of  the  state  deviation  vector  at  time 
tn  is  given  in  equation  (B~8)« 


(B-8) 


The  best  estimate  1  j  the  extrapolated  estimate,  defined  in  equation 
(8-3),  plus  the  weighted  difference  between  the  observed  and  extra¬ 
polated  measurement  deviations*  The  extrapolated  quantity,  d'  de¬ 
fined  as  in  equation  (A-16),  le  what  the  measurement  deviation  is  ex¬ 
pected  to  be* 


(B-9) 


The  vector  ^  io  the  4-dimcnsional  measurement  vector  defined  in 
Appendix  A*  The  weighting  vector,  w  ,  in  equation  (8-8),  io  a  func- 
tion  of  the  covariance  matrix*  The  observed  measurement  c aviation 
differs  fr-xn  ito  trus  value  by  the  meaaurcment  error  a,  a  random 
variable  assumed  to  have  zero  mean  value* 

•  $  «v» ♦  • 


(9-10) 
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equation  (B-3)  is  substituted  into  (B-l)  to  obtain  an  expression 
for  the  error  vector  in  terms  of  known  quantities* 

2*  -  s&'  *4»  ♦*»[*£ 

Sn  -  s&  -  Ss*  ♦$,[*5 *ai>]  <B-“) 

The  vector  £  Is  equivalent  to  S ^  since  the  actual  state  devia¬ 
tion  at  point  n  does  not  change  when  a  new  measurement  is  made* 
Therefore,  in  equation  (B-ll): 


s’ 

“U 


Rewriting  equation  (B-ll): 


e 

—a 


(I  -  Vtt> 


TV  +  vl* 
—a  — <1 


(B-12) 


From  aquation  (B*»2),  the  covariance  matrix  as  a  function  of  ie 
given  in  (8-13)* 

E„  *  (I  •  *2  (B-l3> 

In  equation  (8-13),  the  error  vector  and  the  measurement  error  are 
assumed  independent*  Then  the  cross  terms  a  equal  zero  since  the 
average  measurement  error  is  assumed  zero*  Another  result  of  this 
assumption  in  the  expression  for  ths  variance  of  tho  measurement 
error* 

<T2  *  a2  -  a2  °a2  (Q-14) 

Since  the  weighting  vector,  is  arbitrary,  it  con  be  chosen 
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to  minimize  the  trace  of  the  covariance  raetrix,  thue  minimising  both 
the  position  and  velocity  errors*  Writing  the  perturbation  equation 

i 

of  the  trace  of  (B-13): 

StrfeJ  -  trfiaj  -  tr  [(- £w^)e£(i  .  ^  .  .  . 

+  a  ’  Za*iK(-£aS&  +  <t2  +  SSnSgo-2  -  0  (B-13) 

Since  E*  ia  synrnetric  by  definition,  the  first  two  terms  of  equation 
(B-15)  end  the  last  two  are  the  transpose  of  one  another*  The  fol- 
lowing  relations  hold  for  the  trace: 


tr[ATJ  -  trj/l] 
trj"A  +  b|  •  trj/V  +  BTJ 


Applying  these  relations  to  (B-15): 


2  tr  5  -  &ii£>  ♦  <? 2] 

2  trjs it  [jg  <rz  •  &  <«  *  «£>]!  - 


(B-16) 


Since  the  weighting  vector  is  arbitrary,  the  portion  of  (ml  6)  mul« 
tipi  led  by  $  must  equal  zero* 


«r2  + 


(B-17) 


The  quantity  j^E^  is  a  quadratic  form,  which  makes  it  a  scalar. 
Therefore  the  bracketed  quantity  in  (J-17)  is  a  scalar*  Ihc  trans¬ 
pose  of  the  optimum  weighting  vector  then  has  the  form: 


T 

vL 

-n 


£n  <4 


<r2  + 


(B-18) 
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Since  the  matrix  E*  is  symmetries 


v 
— n 


(B-19) 


The  recursion  formula  for  the  covariance  matrix  is  obtained  by  sub* 
stltutlng  (0*18)  and  (B-19)  back  into  equation  < B— 1 3 > • 

-  [i  .  En&n£n  \  E*  1 1  .  2 

'  ’  T-*  1  r2  +  2nza2m )  ^♦fiXSa)2 


n&i/ 


n 


i^EnSn 

-  “nboiaiEn  „  t-rwnfinan  +  EnSn^EnSn^En 
r2  +  Snt'tfin  <r2  *  ( <T2  +  fi^^,)2 


•  •  • 


P*  n  nTp*  2 

♦ _ :otv3\:p. _ 


En  *  2 


En&nfinEa 


r2+ 


*  <afcfa  +  r2)(E^^) 


T-« 


(4eX 


«■ 


r 


n 


2  +  ff&n 


(B-20) 


The  recursion  formula  for  the  estimation  error  covariance 
matrix  is  somewhat  inconvenient  for  the  steepest  descent  method  of 
solution*  The  Influence  coefficient  expression  obtained  by  using 
(B-20)  would  be  difficult  to  work  with*  A  less  complicated  expres¬ 
sion,  involving  the  Inverse  of  the  En  matrix,  will  be  proved  valid* 
The  proposed  formula  is  given  in  equation  (B*21)* 

E-l  .  «£>-!  * 

<r 


(B-21) 
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if  equation.  (8*21)  if  true,  then  its  product  with  (3-20)  should  equal 
the  Identify  matrix* 


<4  2 


E*  -  E"anjin“n 
n  _  ~  g.2 


VA 


<EVl*&Sa 

<r2 


n 


h  1  h 


h  i 


^EX  +  *- 
2„« 


^  -  (&XK&X 


cr2(&X  ♦  <r2) 


Equation  (B-21)  will  be  used  instead  of  (B-20)  because  of  its  simpler 
form* 


The  estimation  error  covariance  matrix  at  the  target  is  needed 
for  the  cost  function*  Equation  (B-5)  is  used  to  refer  (B-21)  to 
target  co-ordinatca* 


E«  Cen 


(B-5) 


(E.8)'1 


°1 


From  (B-21)* 

(En«)’1 


(C  T)"*p  pT  C 

(C  Vl(E  Vlc  "l  ♦  -K—ffi  > 


an 


n 


on 


<r  2 


-l 


(E  «)  A  -  (E  a)  1  ♦ 


.  i  f  f 
•ax-1  .  -n-n 


n 


n 


T 


(B-22) 


The  vector  f^  io  defined  as  the  4~dimensionel  measurement  vector  in 
target  co-ordinatea*  Equation  (B-7)  provides  a  further  simplifies- 
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tioa« 


If  there  ore  N  measurements  to  be  made*  equation  (B«23)  can 
plied  N  times#  where  Eq  is  the  initial  covariance  matrix* 


<e  ‘r1  ♦ 

o 


From  equation  (B-7)#  oince  the  Nth  moaaurcment  is  the  last! 


<E^'Vl 


(E  °) 

'  N  9 


av"! 


(E0«) 


av-l 


N 


44 


♦  £ 

Urn  1  V 


Henceforth# 
either  E£** 


the  left  hand  member  of  (B-25)  will  bo  referred 
or  A* 


A 


N 

V 

i— 

k«l 


(B-23) 
be  «p- 

(B-24) 

(B-25) 

to  as 

(B-26) 


APPENDIX  C 


STATE  TRANSITION  MATRIX 

In  Stern's  thesis  (Reference  8),  a  general  formulation  of  the 
state  transition  matrix  for  a  two-body  conic  was  derived  end  is  re- 
printed  in  equation  (C«l)*  The  development  was  carried  out  in  the 
pqz  or  flight  path  co-ordinate  system  because  of  the  relative  simplic¬ 
ity  of  matrix  operations*  The  pqr  system  is  a  member  of  t  class  known 
as  reference  trajectory  co-ordinate  systems,  since  its  fundamental 
plane  io  the  plana  of  the  nominal  orbit*  The  system  rotates  about 
its  z*>ax1s,  which  is  defined  perpendicular  to  the  fundamental  plane* 
The  angular  velocity  of  this  rotation  equals  that  of  the  vehicle* s 
velocity  vector*  The  q-axls  lies  parallel  to  the  velocity  vector  and 
the  p-axis  forms  the  right-handed  triad.  Since,  in  a  circular  orbit, 
the  velocity  vector  is  always  perpendicular  to  the  position  vector, 
the  p-axie  lies  along  the  pooition  vector*  Thus,  it  is  obvious  that 
the  Jilight^path  co-ordinate  oystem  is  identical  to  the  more  familiar 
local-vertical  system*  Thia  is  a  fortunate  rcault  since  the  meaoure- 
ment  vectors  developed  in  Appendix  A  ore  more  easily  defined  in  the 
local  vertical* 

For  the  circular  orbit  case,  e  ■  0  in  Stern* s  formula  end  the 
eccentric  and  true  anomalies  are  identical*  In  addition,  the  third 
and  sixth  rows  and  colums  are  eliminated  since  only  the  in-plane 
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8b 


C 


Ji 


MJ1 

NJi 

SJ1 

TJ1J 

MJ1 

-  (c-lb) 

NJ! 

—  (c-lc) 

SJi 

-  (c-ld) 

—  (c-le) 

Symbol 8: 


(c.la) 


e 


n 


eccentricity 
eccentric  anomaly 

7  (EJ  -  El> 

|  (Ej  +  E1) 

mean  angular  motion 
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deviations  from  position  are  to  be  considered*  The  result  is  the 
4X4  matrix  in  equation  (02),  where  fM  ■  fc(fj  -  f^)*  The  factor 
n  is  the  mean  angular  motion,  defined  as  2  divided  by  the  orbital 
period*  Further  simplification  is  possible  since: 

2  sin  £LL-lL  cos  fiJ  -  sln(f j  -  fj)  (C-3) 


Also,  since  the  matrix  will  be  used  to  refer  all  measurements 
to  target  coordinates*  the  J  can  be  replaced  by  a*  Equation  (C-4) 
is  the  resulting  4X4  matrix  . 

The  derivative  of  the  state  transition  matrix  with  respect  to 
time  is  ncaded  in  the  determination  of  influence  coefficients*  This 
derivative  will  be  denoted  Dfll  and  is  given  in  equation  (C-5).  Time 
is  explicit  in  this  formulation  since  nt  ®  f  • 

A  useful  property  of  the  state-transition  matrix  will  be  uti¬ 
lized  to  simplify  computational  procedures*  As  shown  by  Battin 
(Reference  1),  the  state- transition  matrix  is  symplectic,  that  is* 
its  inverse  can  be  computed  by  the  simple  formula: 


C"1  -  -  JCTJ 


/  o  I 

where  J  « 

^•1  0 


(06) 


If  the  C  matrix  is  partitioned  and  the  rule  applied*  the  inverse  of 
C  is  shown  to  be  a  rearrangement  of  its  elements* 


<r 


i 
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<^0 
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The  transpose  of  the  state-transition  matrix  is,  of  course,  also 
symplectic*  In  the  computation  of  influence  coefficients,  the  de¬ 
rivatives  of  both  the  inverse  and  the  inverse  transpose  arc  required* 
The  fact  that  C  and  CT  are  symplectic  greatly  simplifies  these  com¬ 
putations* 

In  the  program,  the  state- transition  matrix  C,  and  its  deriva¬ 
tive  D,  are  computed  according  to  equations  (C-4)  and  (05)  respec¬ 
tively.  The  inverse  and  Inverse  transpose  contain  the  same  elements 
as  C  in  different  positions*  Therefore,  their  derivatives  will  con¬ 
tain  the  same  elements  ac  D  in  thesa  positions  also*  The  easiest  way 
to  form  these  derivatives  is  by  applying  equation  (C-7)  to  the  ele¬ 
ments  of  D  end  Dt*  It  is  important  to  no  to  that  the  results  are  not 
iT*  and  (D*)  but  rather  the  derivatives  of  C"*  and  (CT)***»  The 
formula  for  the  Inverse  of  a  symplectic  matrix  does  not  hold  whsn 
applied  to  D  and  DA,  which  are  not  symplectic*  Equation  (C-7)  is 
merely  tho  easiest  programming  method  which  results  in  the  desired 

-i-CC*1)  and  .l.CC1)*1. 
dt  dt 


quantities: 


APPENDIX  D 


INFLUENCE  COEFFICIENTS  -  TIME  CHANGE 


The  relation  of  the  cost  function  to  the  independent  variable 
is  developed  in  Chapter  3*  For  this  problem,  the  times  of  the  deci* 
•ion  points  ore  varied  and  the  relation  has  the  following  form 

§COBt  -  Z.  £tk  -  S  Lk  Jtk  (D»l) 

k  6  k  k 

The  quantity  Lj.,  the  influence  coefficient,  determines  the  change 
in  cost  due  to  a  change  at  the  k decision  point*  The  relation  is 
accurate  to  the  first  order,  assuming  that  the  time  change  is  small 
enough  so  that  linearity  is  guaranteed*  As  explained  in  Chapter  2, 
the  cost  function  is  the  mean  "squared  arror  in  position  estimation 
at  the  target* 


The  correlation  matrix  of  estimation  errors,  E  »  is  promultlplied  by 
Q  so  that  only  its  first  two  diagonal  elements  are  used  to  determine 
the  cost  function*  For  a  different  cost,  Q  would  have  a  different 
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forest  or  might  not  be  necessary*  It  is  used  only  as  a  convenient 
means  for  expressing  the  cost  function  in  equation  form  and  is  not 
essential  to  the  subsequent  derivation* 

From  equations  (D-l)  and  (D-2),  the  expression  for  the 
influence  coefficient  is  developed  as  follows! 


(D-3) 


An  expression  for  the  error  correlation 
in  Appendix  0* 


N 


matrix  was  determined 


-1 

E  •  A  ■ 

a 


(B-24) 


LQ  ie  the  error  matrix  at  target  due  to  initial  estimation  errors, 

2 

q-  ie  the  varianco  of  the  measurement  noise  and  f^  is  the  modified 
rneosuremant  vector  expressed  in  target  co-ordinates*  Equation  (0-24) 
ia  denoted  A  to  simplify  the  mathenattcc* 


I 


The  partial  derivative  of  EQ  can  be  written  as  follows! 


ijk 


A  ♦ 


CD-4) 


0 
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^  £k 


A" 


1 


if? 

<5  *ie' 


D 


(D-5) 


the  lost  step  following  Lrom  (D~4,* 

Substituting  equation  (D-5)  into  (D-3) x 


«tr 


(D-6) 


An  expression  for  tha  partial  derivative  of  A  is  needed  to  complete 
the  derivation*  Since  the  measurement  vectors  ore  dafined  in  the 
local  vertical  co-ordinate  system  for  each  of  the  decision  points* 
the  state- transit  ion  matrix  must  be  used  to  transfer  them  to  target 
co-ordinate::  This  vas  done  in  Appendix  B,  but  will  be  considered 

here  in  a  slightly  different  manner*  As  noted  in  Appendix  A*  the 
type  of  measurement  used  determines  n  unique  &  vector*  the  measure¬ 
ment  vector  in  local  vertical  co-ordinates  with  two  added  zeroes* 
The  star  elevation  angle  will  be  used  at  each  point,  therefore  the 
k**1  measurement  vector  in  the  two  co-ordinate  systems  are  related  by 
equation  (D- 7)* 

»  -  -  ££  (d-v 

M  le  the  measured  angle  and  x  is  the  state  vector,  where  the  super¬ 
script  a  refers  to  the  target  co-ordlru*te  system*  The  state-transi¬ 
tion  matrix  is  defined  by  equation  (0-8)* 


<  -  c.k4 


(D-8) 


Rewriting  (D-7)j 

«k  2k  "  £  C*K  2k 
"  ik  cak 
£k  “  ^Ak  Ik 

4  - 
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(D-9) 


Each  of  the  N  identical  measurement  vectors  art  transferred  to  tar* 
get  co-ordinates  by  equation  (D-9)* 

Referring  to  equation  (B*5): 


(B-5) 


Applying  (B-5)  to  the  Initial  error  matrix,  E®,  equation  (B-24)  can 
now  be  written  as  an  explicit  function  of  time* 


A 


<w£>*  ♦  1 

k®i  <r  2 


(D-10) 


The  influence  coefficient,  equation  (D-6),  requires  the  partial  deriv¬ 
ative  of  A  with  recpcct  to  tjc* 


^  A 

£  4 


p(cj)-1 

L 


£|<£kCak  +  (C»k) 


(D-ll) 


A  simplified  expression  for  the  influence  coefficient  can  now  be  ob¬ 
tained,  using  the  following  identities  involving  the  trace* 

tr  [ab]  o  tr  [baJ 
tr  [a  ♦  b]  -  tr[A  ■*■  BTJ 
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From  (D-6)j 


Lk  ■  *' 


h  it  E-1 


-  -trri* 

L  an 


"  r2  tr  [(~(C'V"'  +  —f)  V*.  ] 


p  tr  [-- &*&  c*  E*QE«] 


p tr  fc  cIk  Ee^a  J 


(D-12) 


The  expression  In  brockets  In  equation  (D«12)  is  in  the  quadratic 
form,  hence  it  is  a  scalar  quantity  and  the  trace  notation  can  b« 
dropped*  The  final  expression  for  the  kta  influence  coefficient  is 


given  in  (I>-13)* 


T  ck 


(D-13) 


APPENDIX  E 


INFLUENCE  COEFFICIENTS  -  HOIUZON  CHANGE 

The  Influence  coefficients  associated  vith  a  possible  change 
of  horizon  have  the  same  function  ea  those  sssoclsted  vith  the  time 
change  and  the  cost  relation  has  a  similar  form: 

£  coit  -  Z  § 4  ■  ?  i,,  ^  (e-1) 

k  k 

is  defined  in  equation  (E-2),  where  ^  ia  the  measure- 

th 

meat  vector  expressed  in  the  local  vertical  syatem  at  the  k  point* 

i^(HEW)  -  ^(OLD)  +^"4  <E-2) 

The  Influence  coefficient  L^,  in  equation  (*-l)»  la  a  column  vector 
whose  scalar  product  vith^  is  the  cost  change  resulting  from  a 
change  in  horizon  at  the  k^  point*  As  indicated  in  (E«l),  the  sum 
of  tha  k  cont  changes  should  result  in  the  total  cost  change,  if  the 
assumption  of  linearity  is  valid*  The  some  cost  function,  position 
estimation  error,  is  used  and  so  the  Q  matrix,  whose  purpose  was 
explained  in  Appendix  D,  is  retained* 

Cost  -  tr  [QE^  (E-3) 

th 

The  expression  for  the  cost  change  at  the  ktn  point  is  found  by 


The  vector^* 
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writing  the  perturbation  equation  of  <e-3)  for  a  change  of  horiaon* 


■  f. . 

<  &  Coat>k  -  £  t 

(  £  Cost),  •  tr  < 
k 

k4 

(  §  Cost)k  -  tr 

m 

(E-4) 


Equation  (D"10)  provides  an  egression  in  terms  of  known  quantities 
for  the  inverse  of  the  estimation  error  covariance  matrix* 


.-1 

"a 


A  ■ 


<CaoEoCJ>*1 


+ 


N 

I 

k»l 


<C 


Tj-l 


T„  -1 


Mh.flk, 


(D-10) 


The  perturbation  of  Efi  can  be  written  as  a  function  of  A*  similar  to 
the  development  in  Appendix  D* 

...  v  -  I 

$  EaA  +  sa  s  A  •  0 

S \  m  -  A*1 

S'5,  ■  *  Ea^  A®,  (E-J) 

•farrlnfi  to  CE-4)  ,  the  it***  cost  change  now  has  the  form 

Co»t)k  -  .  tr  [qe.  £  A  E.J  (E*6) 

The  trace  of  a  product  of  matrices  can  be  rearranged  aa  follows! 

trfABj  -  tr[k]  (e-7) 

.  ’  •  .  •  t  ‘  .V  .  ■  ’  Z  / 

Applying  (E«7)  to  <E-6)»  with  "A"  in  (E-7)  equal  to  E  ,  results  in 

A 
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more  convenient  expression* 


(  g  Coat)k  •  -  tr  Je^QE^aJ 


(B*8) 


the  perturbation  of  A»  at  the  kth  point,  can  be  written  directly 

equation  (D-10). 


k  [(CJ>*V  Mk  c^1]  CE-»> 


rr?- 


thc  kw  cost  change  can  now  be  written  as  a  function  of  known  quan- 


•it ins  from  equations  (E-8)  and  (E»9)« 


$  Cost)k  «  -  ^2  fcr  i  Ea^Ea 


[(cak>  ^£kfik  Cakl  +  (Cak*  ^ 


(E-10) 


Separating  the  trace  expression  in  (E-10)  into,  two  parts  ^nd  applying 


(E«7)  with  ’’A"-  equal  to  EaQEfl: 


(  o  Cost)k  -  .-lj  iM^apsk* CakX 


tr  Si  C^1]  EaQE 


(E-ll) 


.toother  identity  involving  the  trace  is  found  to  be  useful* 


tr  fV]  -  tr  J^aJ 


.  CE-12) 


vn  examination  of  aquation  (E-ll)  reveals  the  following  identltyi 


v5®*]*  -  (e.13) 


Sines  the  EQ  and  Q  matrices  are  both  ayrrroetric,  the  "T"*e  on  these 
matrices  can  be  dropped*  Equations  (E-12)  »*nd  (E-13)  are  then  applied 
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to  to  further  reduce  the  cost  change  expression* 


(  S  dost), 


**  -  JL  tr  ,/E 

<r2 


««Ea[<CJ 


Tv**  T  ^  *l] 

5  44  c«k  J 


(B-14) 


The  elerasnto  of  the  trace  can  be  rearranged  by  twice  applying  equa* 
tlon  (E-7)  vhsre  "A"  in  (E-7)  is  first  equal  to  E  QE  and  then 


<  $  Cost), 


-2  tr  <&£  Cak 


lEaQE«(cak)*1  fikj 


(£•15) 


Equation  (E**15)  can  be  simplified  by  considering  ths  following  group* 
inga  of  its  elemental 


c  ?}  e  qe  (c  £)*a 
ck  ay  av  ak' 


s  4  X  4  matrix 


a  row  vector 


Since  ^  is  a  colum  vector,  the  bracketed  expression  in  (E-15) 
is  o  scalar  and  the  trace  notation  can  be  dropped* 

(  <5  05oOk  -  .-Zj  C^1  E»QE4(cJ)*V  £^j  (E-16) 

The  influence  vector  may  be  identified  by  comparing  equation  (E-16) 

vith  (E-2).  .  . 


C  '  -^l[4TcAlEaQE<.(c«k)*lJ 


(E-17) 


In  the  time  change  portion  of  the  problem,  the  time  of  the 
decision  points  could  be  changed  In  whichever  direction  resulted  in 
ffl  decrease  in  cost*  There  Is  no  such  control  in  changing  the  hori¬ 
zon  sines  th®£  ^  vector  is  predetermined  by  the  present  horison* 
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There!  arc  only  evo  possible  value*  for  $  depending  on  whether 
the  proposed  change  Is  from  the  left  to  right  horison  or  vice- verse* 


'rom  Appendix  At 


%  (LEFT)  » 


4 


(RIGHT) 


a (a2  -  r  2)^ 

E 

1 

» 

S 

0 

0 


(A-H) 


CA-14) 


The  new  measurement  vector  is  related  to  the  old  by  equation  (E-3). 
Therefore s  the  t vo  £"  vectors  Hava  the  following  form* 


S’  (LEFT  TO  RIGHT)  ■= 


(E»18) 


g  g^  (RIGHT  TO  LEFT) 


(E*19) 


0 
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As  explains  la  Chapter  2,  equation  (E-20)  is  applied  at  each  point 
and.  If  the  resultant  change  in  coat  is  negative,  the  proposed  hori¬ 
zon  change  is  made* 


(  S  Oost)^  - 


(E-20) 


APPENDIX  F 


COMPUTER  PROGRAM 


The  Fortran  program  used  to  teat  the  theory  covered  la  this 
thesis  is  quite  complex  and  involves  many  subroutines*  The  main 
program  flow  charts  ore  shown  in  Chapter  4*  Also  the  logic  for 
three  subroutines  is  covered  in  Chapter  4*  The  main  program  carries 
the  primary  load  for  computing  the  new  measurement  time  schedule, 
selecting  the  measurement  horiron  schedule,  and  performing  the  opti* 
mi  nation*  The  subroutines  compute  various  matrices,  perform  various 
matrix  operations,  create  new  measurement  time  and  ho risen  schedules, 
and  print  various  matrices* 

The  input  data  needed  by  the  main  program  has  the  following 

sequence  and  unltsi 

1)  Variance  of  the  measurement  (radians^) 

2)  Altitude  of  spacecraft  (miles) 

3)  Radius  of  the  earth  (mi lea) 

4)  Final  angle  (degrees) 

5)  Gravitational,  constant  (feet^/sec^) 

6)  Kumber  of  iterations 

7)  Initial  measurement  time  schedule  (seconds) 
so  Block  one,  right  then  left 
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b«  Stock  two,  initial  schedule 

8)  Initial  estimation  error  covariar.ee  matrix  (feet2  end 
£®>2t2/ second2) 

9)  Maximum  time  changes  (seconds) 
a*  Right  then  left 

10)  initial  measurement  horizon  schedule  (Block  two  only, 

right  horizon— flag  equals  one,  left  horl son— flag  equals 

sere) 

The  standard  matrix  operations  are  needed  in  the  program*  Also 
since  sem®  of  the  vectors  and  matrices  depend  upon  the  particular 
m^aRirsnent  time,  special  variations  are  needed*  The  operations 
needed ,  and  the  subroutines  performing  them  ares 

1)  Dot  product  (MULTF  and  MULTG), 

2)  Dlatic  product  (MULTA  and  MULTAA), 

3)  Matrix  inversion  (INVERT), 

O  Matrix  multiplication  (MULTB,  MULTC,  MULTD,  and  MULTE), 

5)  Product  of  row  vector  and  matrix  (MULTI), 

6)  Matrix  sums  (SUMA,  SUMB,  and  SUMC), 

7)  Trace  (TRACS  and  TRACEB)* 

Various  matrices,  vectors,  and  the  results  of  several  matrix 
and  vector  equations  are  needed  in  the  program* _ These  arcs 

1)  State  transition  matrix  (TRANSK), 

2)  Inverse  of  state  transition  matrix  (TRANMI), 

3)  Transpose  inverse  of  state  transition  matrix  (TRANIT), 

4)  Derivative  or  each  element  of  the  state  transition  matrix 
(DTRANM), 


109 


5)  Matrix  Q  (AQM), 

6)  Initial  estimation  error  covariance  matrix  in  target 
ordinates  equation (  B*3  )  (ERRORI), 

7)  Target  estimation  error  oo variance  matrix,  aquation 

(B-26)  (ERRORM), 


8)  Horison  vector  matrix  (HORIZS), 

9)  Horizon  change  vector,  (LAMCA), 

10)  Equation  F»1  (DERRGR), 


Tv  "1 


(P-1) 


j  <cak> 

^  b. 


T  -  -1 


cak 


er¬ 


ic  o  T  ^  ^Cak^"  ~ 
'CaJr  o\ 

2 


11)  Teat  on  time  change  (CHECK), 

12)  Hm  measurement  time  schedule  (LOGIC), 

13)  Nm*  measurement  horieon  schedule  (JUMP). 


Various  matrices,  vectors,  lists,  and  headings  are  needed  in 

tha  results.  These  are  printed  by  PRLNTA,  PRINTB,  PRINTC,  PR  INTO, 

and  PRINTR* 


Printouts  will  follow: 


r\  n 


MAIN  BLOCK  ONE 


f  L  1ST 

^  LABEL 

definition  of  the  variables  and  constants 

w  is  the  angu  ar  velocity  of  the  spacecraft 

to  is  the  TIME  of  orbit 

fa  is  the  final  angle  or  target  angle 

T(I)  is  the  initial  timing  schedule 

El  IS  THE  initial  ESTIMATION  ERROR  COVARIENCE  MATRIX 

DEL Tmr  AND  DELTML  ARE  THE  MAXIMUM  TIME  CHANGES 

sigmas  is  the  variance  of  the  measurement 
al  is  the  altitude  of  the  craft 

RE  IS  the  RADIUS  OF  THE  EARTH 

U  JS  THE  GRAVITATIONAL  CONSTANT 

HR  AND  HL  ARE  THE  MEASURMENT  VECTORS 

EAR  AND  EAL  ARE  THE  FINAL  CORRELATION  MATRICES 

DIMENSION  DR( 4.4) »DL<4»4)  »DELTR(30)  »DELTL  <  30  )  *  EAL  <  4*4  )  * 

1  EAR(4*4) *EEL( 30*4*4) »EERf30*4*4)»EI(4*4)» ERRL  ( 30  *4*  4 ) * 

1  ERRR( 30*4*4)  *  HHL (4*4) »HHL(4»4) »HL ( 4 ) *HR ( 4 ) *Q ( 4 *4 ) , 

3  RR(30,4*4  )  *SS( 30*4*4) »Tal( 30 ) .TALI ( 30) .TAR (30) »TARI (30 ) . 

4  T ( 3 0 )  *  TR ( 3 0 )  »  TL ( 30 ) 

common  fa* w» sigmas *t*ei .deltmr* deltml 

1000  read  1 »SIGMAS*AL*RE»FA 
I  FORMAT ( F20 . 10 ) 

111  READ  1 1 2  »  U 

112  FORMAT (E20. 10} 

PRINT  800* SIGMAS 

300  FORMAT  ( 1H0 1 10X34H  THE  VARIANCE  OF  THE  MEASUREMENT  «.F20.10 
PRINT  8  0 1  *  AL  trcv.iu 

B°1  FORMAT ( 1  HO  *  1 OX 3 3H  THE  ALTITUDE  OF  THE  SPACECRAFT  ®»F20.10. 
16H  MILES) 

PRINT  802.  RE 

802  FORMAT (1H0.J0X26H  THE  RADIUS  OF  THE  EARTH  -.F20.10.6H  MILES 

PRINT  8  0  3  ♦  F  A  — 

303  E0RMAT(1H0,10X18H  THE  FINAL  ANGLE  -.F20.10.8H  DEGREES) 

PRINT  804.U  ’ 

804  FORMAT! 1H0,10X31H  THE  GRAVITATIONAL  CONSTANT  U  =*E20.10* 

1 3 OH  FEET  CUBED.  PER  SECOND  SQUARED) 

READ  999 *NT 
999  FORMAT  (  HO) 

PRINT  998, NI 

998  FORMAT ( 1H0, 10X28H  THE  NUMBER  OF  ITERATIONS  IS.I6) 

READ  2 » ( TR ( I ) *  1  =  1 .30 ) 

READ  2. ( TL( I )  .  I»1  .30 ) 

PRINT  808 

803  FORMAT ( 1H5 » 10X36K  THE  INITIAL  TIMING  SCHEDULE  FOR  THE. 

1  17H  RIGHT  HORIZON  IS) 

CALL  PRINTA(TR) 

PR  T  NT  388 

888  FORMAT ( 1H5 * 10X36H  THE  INITIAL  TIMING  SCHEDULE  FOR  THE. 

1  16H  LEFT  HORIZON  IS) 


NUMBER  OF  ITERATIONS  IS.I6) 


INITIAL  TIMING  SCHEDULE  FOR  The. 


INITIAL  timing  schedule  for  the. 


Ill 


2 

2000 


CALL  PRINTA(TL) 

FORMAT  (F20.10) 

READ  2000# { (El (j#<) # j m i # 4 ) 
FORMAT  (F25.8) 

PRINT  807 


807  F0RMAT(1H5,10X28H  The  INITIAL  ERROR  MATRIX  IS) 

CALL  PRINTB  (El) 

READ  2001,  DELTMR #DELTML 

2001  FORMAT  (F20.10) 

PRINT  2100, DELTMR 

2100  FORMAT! 1H0.10X29H  TIME  INCREMENT  RIGHT  HOR I ZON . F20. 10 ) 
PRINT  2101, DELTML 

2101  FORMAT ( 1H0 , 10X28H  TIME  INCREMENT  LEFT  HOR I ZON # F20 . 10 ) 

AL  =  AL*5  2  80 , 


RE=RE*5280. 

F  A  =  F  A*  3  *  1 4 1 5 9  / 1 8  0  , 

3  T 0=2, 0*3,1 4159# ( ( RE+AL ) **1 . 5 ) /SQR T F ( U ) 

4  W  =  2 • * 3 • 14159/TO 
PDINT  805, TO 

805  FORMAT ( 1H0 » 1 0X20H  THE  TIME  OF  ORBIT  =»E20«10»8H  SECONDS) 
PRINT  806, W 

806  FORMAT ( 1H0 , 10X41H  THE  ANGULAR  VELOCITY  OF  THE  SPACECRAFT  *. 
1E20.10, 1 9H  RADIANS  PER  SECOND) 

5  HR(1 )=RE/( ( AL+RE)*SQRTF( ( ( AL+RE ) **2 ) - ( RE**2 ) ) ) 

6  HR(2)=-1*/ (AL+RE ) 

7  HR ( 3 ) =0  « 

8  HR ' 4 ) =0 , 

9  HL ( 1 ) =  +  HR ( 1 ) 

10  HL ( 2 ) =-HR ( 2 ) 

1 1  HL (3 )»0  * 

12  H  L  (  4 ) *  0  , 

PRINT  809 

809  FORMAT ( 1H5 #10X28H  THE  RIGHT  HORIZON  VECTOR  IS) 

PRINT  8  8  0#  HR 

PRINT  810 

810  FORMAT ( 1H5 » 10X27H  THE  LEFT  HORIZON  VECTOR  IS) 

PRINT  880,  HL 

880  FORMAT  ( 1H0# 4E30.10) 

N  I  I  =0 , 0 

13  CALL  ERRORM { EAR ♦ HR » OLCOSR  » TR } 

14  CALL  ERRORMf EAL » HL , OLCOSL t TL > 

PRINT  811 

811  FORMAT(1H5,10X40H  FINAL  CORRELATION  MATRIX#  RIGHT  HORIZON) 

CALL  PRINTB  (EAR)  '  "  ' 

PRINT  812 

812  FORMAT( 1H5,10X39H  FINAL  CORRELATION  MATRIX#  LEFT  HORIZON) 

CALL  PRINTB  (EAL) 

T  A  R  I  M  =  0  „  0 
TALIM  =  0 ,0 

15  DO  341  1=1,30,1 

16  CALL  AOM(O) 
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17  CALL  MULTB  (DR.EAR.Q) 

18  CALL  MULTB  (  DL *  EAL  *Q ) 

19  CALL  DERROR { ERRR  »HR ♦ I *TR ) 

21  CALL  DERROR(ERRL*HL* I »TL ) 

CALL  MULTF  (  RR , £RRR , DR  .  I  ) 

210  CALL  MULTE  ( SS , ERRL , DL . I) 

22  CALL  M U L  T C  ( EER . EAR , RR , I ) 

220  CALL  M'JLTC  (  EEL  »  EAL .  SS  .  I  ) 

23  CALL  TRACEB  (  TAR  ♦ EER  t I  ) 

24  CALL  TRACEB  (TALtEELtI) 

25  TAR I( I  5  =-T  AR ( T  )*W 

26  TALI ( I)»-TAL( I)*W 

29  Thru  34  DETERMINE  the  LARGEST  INFLUENCE  COEFFICIENT 
FOR  both  the  left  and  right  horizon 

TARIMtRIGHT  AND  TALIM.LEFT 

29  I F ( ABSF ( TAR  I <  I ) ) -ABS F ( T AR I M ) ) 311 *33*33 
311  GO  TO  30 

33  T  ARIM«TARI  ( I  ) 

30  IF (ABSF (TALI ( I) ) -ABSF ( TALI M ) 5321.34,34 
321  GO  TO  341 

34  T  AL IM  =  T  AL I ( I ) 

341  CONTINUE 

PRINT  261 

261  FORMAT ( 1H0 » 1OX37H  RIGHT  HORIZON  INFLUENCE  COEFFICIENTS) 

DO  264  1  =  1  ,30,1 

26  -  PRINT  263 , TAR  I  ( I  )  »  I  “  •  — 

PRINT  262 

262  FORMAT ( 1H0 , 10X36H  LEFT  HORIZON  INFLUENCE  COEFFICIENTS) 

DO  265  1=1,30,1 

265  PRINT  263, TALI (  I  )  ,  I 

263  F0RMAT(1H0,10X12H  COEFF I C I EN T t E2 0 . 1 0 ♦ 5X 1 2H  TIME  NUMBER. 14) 

36  PRINT  362 

362  f0RMAT( 1H5.10X23H  DATA  FOR  RIGHT  HORIZON) 

361  CRUD=DELTMR 

35  CALL  CHECK(TARIM»TARI,OLCOSR»HR*TR,CRUD) 

37  CALL  LOGIC(BIGA»TARIM*CRUD,TARI » OLCOSR » PCOSr , EAR » HR  * TR ) 

38  PRINT  382 

382  FORMAT ( 1H5 » 10X22H  DATA  FOR  LEFT  HORIZON) 

381  CRUD=DELTML 

CALL  CHECK  ( T AL IM . T AL I .OLCOSL * HL * TL * CRUD ) 

39  CALL  LOGIC(BIGB»TALIM,CRUD*TALI .OLCOSL, PCOSL, EAL, HL*TL1 

40NJI*NII+1 

41  IE(NII-NI)  42.45.45 

42  PRINT  43,  Nil 

^3  FORMAT ( 1H0 , 10X15H  THE  END  OF  THE.I5.10H  ITERATION) 

44  GO  TO  13 

45  PRINT  46 

46  FORMAT ( 1H0,10X8H  THE  END) 

47  GO  TO  1000 
END 


n  n 
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main  block  two 


LIST 


LABEL 


DEFINITION  OF  THE  VARIABLES  AND  CONSTANTS 
TO  IS  THE  TIME  OF  ORBIT 

W  IS  THE  ANGU  AR  VELOCITY  OF  THE  SPACECRAFT 
FA  IS  THE  FINAL  ANGLE  OR  TARGET  ANGLE 
T(I)  IS  THE  INITIAL  TIMING  SCHEDULE 

El  IS  THE  INITIAL  ESTIMATION  ERROR  COVARIENCE  MATRIX 
DELTMR  and  deltml  ARE  the  maximum  time  CHANGES 
sigmas  is  the  variance  of  the  measurement 
AL  is  THE  altitude  OF  THE  CRAFT 
ihv(i)  is  the  initial  horizon  schedule 

RE  IS  THE  RADIUS  OF  THE  EARTH 
U  IS  THE  GRAVITATIONAL  CONSTANT 
HR  AND  HL  ARE  THE  MEASUREMENT  VECTORS 

ea  is  the  final  correlation  Matrix 

DIMENSION  DELTR(30),DHLR(4)  » DHRL ( 4 ) • FA ( 4 * 4 ) * E I < 4 * 4 ) » 

1  EER ( 30  »4»  4 ) *  ERRR ( 30  .4.4 ) *HR ( 4 )  * HL ( 4 ) * IHV ( 30 ) * 

2  HS  {  30  » 4  )  »DELCOS(30)  *Q  (  4.4)  »RR(30»4»4)  »TAR_(30)  *  T ARI { 30 )  * 

3  TLAMDA ( 30 ) ,TE30) 

COMMON  FA»W»SIGMAS»T  »EI r. DELTMR * DELTML 
1000  READ  l.SIGMAS.AL.RE.FA 
1  FORMAT ( F20, 10  ) 

111  READ  1 1 2  t U 

112  FORMAT (E20. 10) 

PRINT  8  00  *  S I GMAS 

800  FORMAT  ( 1H0 ♦ 1 0X34H  THE  VARIANCE  OF  THE  MEASUREMENT  -.F20.10) 
PRINT  801 »al 

801  FORMAT  I 1«0,10X33H  THE  ALTITUDE  OF  THE  SPACECRAFT  -.F20.10, 
16H  MILES) 

PRINT  80  2 ♦  RE 

802  FORMATE 1H0.10X26H  THE  RADIUS  OF  THE  EARTH  -iF20.10.6H  MILES) 

PRINT  803 ♦ FA 


8O3  r ORMAT ( 1H0 ♦ 1 0X1 8H  THE  FINAL  ANGLE  -.F20.10.8H  DEGREES) 
PRINT  804 ♦ U 


804  FORMAT ( 1H0.10X31H  The  GRAVITATIONAL  CONSTANT  U 
1 3 OH  FEET  CUBED  PER  SECOND  SQUARED) 

READ  999 .NT 
999  FORMAT  (110) 


“.E20.10. 


PRINT  998. NT 

998  FORMATE 1H0.10X28H  THE  NUMBER  OF  ITERATIONS  IS.I6) 
READ  2 ♦  ( T (  I  )  ,  1*1 ,30 ) 

2  FORMAT  (F20.1C) 

PRINT  808 

808  FORMATE 1H5.10X31H  THE  INITIAL  TIMING  SCHEDULE  IS) 
CALL  PRINTA  (T) 

READ  2000,((EICJ»K)»J«1»4)*K«1*4) 

2000  FORMAT  EF25.85 
PRINT  807 
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807  FORMAT!  1H5 . 10X28H  THE  INITIAL  ERROR  MATRIX  IS) 

CALL  PRINTS  (El  ) 

READ  2001*  DELTMR*DELTML 
200T  FORMAT  (  F20.10) 

PRINT  2 1 00  *  DELTMR 

2100  FORMAT  (  1  HO  * 10X15H  TIME  INCREMENT  * E20 . 10 ) 

READ  2002* ( IHV(  I  )* 1*1 .30) 

2002  FORMAK 15) 

PRINT  2003 

2003  FORMAT  (  1 H5  *  1 0X32H  THE  INITIAL  HORIZON  SCHEDULE  IS) 

CALL  PR  I NTD ( IHV) 

AL=AL*5280. 

RE*RE*5280. 

FA*FA*3. 14159/180* 

3  T0*2«0*3 .1 41 59* ( (RE+AL ) *  * 1  •  5 ) /SQRTF(U) 

4  W*2 .*3. 14159/TO 
PRINT  805. TO 

805  FORMAT  (  1  HO  t 1  OX 2 OH  THE  TIME  OF  ORBIT  *.E20.10.8H  SECONDS) 
PRINT  806. W 

806  FORMAT  (  1H0*  10X.41H  The  ANGULAR  VELOCITY  OF  THE  SPACECRAFT  *» 
1E20.10.19H  RADIANS  PER  SECOND) 

5  HR ( 1 ) =RE / ( ( AL+RE)*SORTF( (  (  AL  +  RE) **2 )-(RE**2 ) ) ) 

6  HR( 2  )  *  - 1 •/ ( AL+R  E ) 

7  HR  (  3  )  *0  . 

8  HR ( 4  )  *0 . 

9  HL ( 1 )  *-HR ( 1  ) 

10  HL ( 2 ) =HR ( 2  ) 

11  HL ( 3  )  =0 . 

12  HL ( 4 ) *0 . 

PRINT  809 

809  FORMAT  (1H5.10X28H  THE  RIGHT  HORIZON  VECTOR  IS) 

PRINT  880. HR 

PRINT  810 

810  FORMAT(  1H5.10X27H  THE  LEFT  HORIZON  VECTOR  IS) 

PRINT  880.  HL 

DHRL(1)=2.0*HL(  1  ) 

DHRL ( 2 ) *0*0 
DHRL ( 3  )  *0*0 
DHRL ( 4 ) *0.0 
DHLR ( 1 ) * 2 •  0*HR  (  1 ) 

DHLR(2) *0.0 
DHLR (3) =0.0 
DHLR (4) =0.0 
PRINT  820 

820  FORMAT!  1H5. 10X40H  THE  DELTA  VECTOR  FROM  R  To  L  HORIZON  IS) 
PRINT  880 ♦ DHRL 

PRINT  821 

821  FORMAT!  1H5.10X40H  THE  DELTA  VECTOR  FROM  L  To  R  HORIZON  IS) 
PRINT  880. DHLR 

880  FORMAT!  1H0.4E30.10) 

N  1 1  *0 


13  CALL  HOR IZS( IHV*HR*HL*HS ) 

CALL  ERRORM(  EA*HS*OLCOSR  ) 

PRINT  811 

811  FOPMAT<lH5»lOX25H  FINAL  CORRELATION  MATRIX) 

CALL  PRINTB(EA) 

TARI M  =  0  *0 

15  DO  341  1*1* 30*1 

16  CALL  AOM(O) 

17  CALL  MULTB(D*  EA  *0 ) 

19  CALL  DERROR ( ERRP * HS* I ) 

21  CALL  MULTCtPR .ERRR.D. I) 

22  CALL  M’JLTC(EERtEA.RR  .  I) 

23  CALL  TRACER  (TAR  *  EER • I ) 

25  TARI ( I ) =-TAR (  I ) *W 

C  TARKI)  ARE  THE  INFLUENCE  COEFFICIENTS 

C  29  Thru  341  DETERMINE  the  LARGEST  INFLUENCE  COEFFICIENT 

29  I F(AB5F(TAR1 ( I) ) -ABSF( TARIM) ) 311 *33*33 
311  GO  TO  341 
33  TARIM*TARI U ) 

341  CONTINUE 

DO  263  1*1*30*1 
263  PRINT  26  1  .  TAR  I  (!)  *  I 

261  FORMAT(2HO,lOX23H  INFLUENCE  COEFFICIENT  *E20*l0* 

1  5X12H  TIME  NUMBER. 13) 

35  CALL  CHECK  ( T AR I M .TAP  I .OLCOSR *HS ) 

361  CRUD*DELTMR 

37  CALL  LOGIC  (  B  IGA  *  TAR  I M *CRUD* T AR I  * OLCOSR . PCOSR * EA * HS ) 

38  CALL  JUMP  ( TL AMDA *DHLR *DHRL * PCOSR . IHV .HR . HL ) 

PUNCH  462  *  IHV 

462  F0RMATU5) 

40  N I  I SN I  1  +  1 

41  IF(NII-NI)  42  *45  *  45 

42  PRINT  43. Nil 

43  F0RMAT(1H0,10X15H  THE  END  OF  THE*I5.10H  ITERATION) 

44  GO  TO  13 

45  PRINT  46 

46  FORMAT ( 1H0.10X8H  THE  END) 

47  GO  TO  1000 
END 


o  u  o  o  o  o  o  o  ^  o  ooo  a  a  o  c 
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SUBROUTINE"  INVERT  (N*QQ) 

C  REFERENCE ( 7 ) 

♦  LIST 

LABEL 

subroutine  inverts  matrix  by  simultaneous  DOUBLE  PRECISION 
ROW  REDUCTION  OF  THE  MATRIX  TO  IDEM  AND  IDEM  TO  THE  INVERSE 
DIMENSION  00(4*4 ) *0(4*8 ) 

DO  10  1*1*4 
DO  5  J= 1  *  4 
0(1*  J)=QQ(I»  J) 

0 ( I  * J  +  4)*0.0 
0( I  * J+8)=0.0 
5  0(1*  J+12 ) *0*0 

10  Q( I  *  1+4) *1  .0 
DO  30  1*1*  N 
DO  14  J*I.  N 

I F ( ABSF (0(1*  I ) ) -ABSF (0 ( J*  I)))  11*  14*  14 
TEST  COR  LARGEST  ELEMENT  IN  COLUMN 

11  DO  12  K*1 *  N 
S*Q(J*  K ) 

0(J*  K ) *0 (  I  *  K) 

0(1*  K)»S 
S*Q( J,K+4 ) 

0(J*K+4)*0( I *  K+4  ) 

12  0 ( I  *  K+4 ) *S 

TRANSFER  ROW  OF  largest  ELEMENT  TO  FIRST  ROW 

14  CONTINUE 
D I V*Q ( 1*1) 

DO  15  J* 1  *  N 
0  (  I  *  J ) *0 ( I  ♦  J ) /D I  V 

15  0( I  *  J+4 ) *0 ( I  *  J+4  )  /DI V 
DIVISION  BY  DIAGONAL  ELEMENTS 
DO  30  J* 1  *  N 
IF(I-J)  20,30*20 

20  D I V=Q ( J  *  I ) 

DO  25  K*  1  *  N 

Q(J*K)*Q(J*K)-Q(I*K)*DIV 
25  Q ( J*K+4 )=Q( J*K+4 )-0( I  *  K  +  4 ) *D I  V 
DIAGONALIZATION  OF  MATRIX 
30  CONTINUE 
DO  35  I  *  1  *  N 
DO  35  J  =  1  *  N 
35  00( I t J)*0(  I*J  +  4) 

RETURN 
END 
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SUBROUTINE  MULT  A ( HH  *  H ) 

*  L  1ST 

*  LABEL 

this  subroutine  performs  the  multiplication  of  a  row 
C  VECTOR  by  a  column  vector 

DIMENSION  HH  (4*4)*  HU) 

1111  DO  1 1 1 A  I  *  1  *  4  *  1 

1112  DO  1114  J*  1 .4 •  1 

1113  HH ( I »J  >*0.0 

1114  HH ( I  . J)*H(  I )#H(  J) 

RETURN 

FND 


SUBROUTINE  MULT AA ( HH . Ht  I) 

*  LIST 

*  LABEL 

c  this  subroutine  performs  the  multiplication  of  a  row  vector 

C  BY  A  COLUMN  VECTOR.  EACH  VECTOR  DEPENDS  ON  THE  TIME  I. 

DIMENSION  HH (30*4*4) *H( 30*4  ) 

1  DO  4  J»1 ,4 .1 

2  DO  4  K* 1 » 4  *  1 

3  HH(  I  *J.K)«=0.0 

4  HH ( I »J.K)«H( I  » J ) *H (  I *K) 

RETURN 

END 


SUBROUTINE  MULTB(XX*XY*XZ) 

*  LIST 

*  LABEL 

C  THIS  SUBROUTINE  MULTIPLIES  TWO  MATRICES 

DIMENSION  XX (4*4).  XY(4*4)*  XZ(4.4) 

2111  DO  2115  J* 1  *  4 . 1 

2112  DO  2115  K-l.4.1 

2 1 13  XX(J.K)*0.0 

2114  DO  2115  L* 1 *4 *  1 

2115  X X ( J » K ) =XX(J*K)+XY(J*L)*XZ(L*K) 

RETURN 

END 


u  u  #  *  U  U  O  O  * 
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SUBROUTINE  MUlTC < WX * WY ♦ WZ . I ) 

*  LIST 

*  LABEL 

THIS  SUBROUTINE  MULTIPLIES  TWO  MATRICES 
THE  result  ANO  multiplicand  INVOLVE  THE  I  TH  TIME 

DIMENSION  WX( 30*4*4 ) *WY<  4*4) *WZ ( 30.4, 4) 

3111  DO  3115  J* 1 #  4  #  1 

3112  DO  3115  X* 1 # 4 # 1 

3113  WX (  I  * J  *X ) *0.0 

3114  DO  3115  L* 1  *  4  #  1 

3115  WX( I#J*X)*WX< I *J*X)+WY(J*L)*WZ( I#L*X> 

RETURN 
END 


SUBROUTINE  MULTD ( UX * UY # UZ . I ) 


4111 


4112 

4113 

4114 

4115 


LIST 

LABEL 

THIS  SUBROUTINE  MULTIPLIES  TWO  MATRICES 

The  RESULT.mUlTIPLIER.AND  MULTIPLICAND  INVOLVE  THE  I  TH  TIME 
DIMENSION  UX( 30*4*4) *UY( 30*4*4) #U2( 30*4*4) 

DO  4115  J  *  1 ♦ 4  # 1 
DO  4115  X*1 *4*  1 
UX ( I • J • X ) *0.0 
DO  4115  L  *  1 »  4  *  1 

UX<  I#J*X)*UX(  I  *  J  *  X ) +U Y (  I  *  J  *L  )*UZ ( I »  L  *  X ) 

RETURN 

END 


SUBROUTINE  MULTE  <RX#RY*RZ*I) 


THIS  SUBROUTINE  PERFORMS  THE  MULT ICL ICAT ION  OF  TWO  MATRICES 
THE  RESULT  AND  MULTIPLIER  INVOLVE  THE  I  TH  TIME 
LIST 
*  LABEL 

DIMENSION  RX( 30*4*4} *RY( 30*4*4) *RZ (4*4) 

561  DO  565  J  =  1 • 4 *  1 

562  DO  565  X=l*4*l 

563  R  X  (  I  *  J  *  X  )  * 0  . 

564  DO  565  L  =  1  *4* 1 

565  R  X ( I  * J  *X  )  *RX (  I  * J»X)+RY<  '  * J * L > *RZ < L *  X ) 

RETURN 

END 


SUBROUTINE  MUlTF(A»B*C) 


LIST 

LABEL 

this  subroutine  performs  the  dot  product  of  a  vector 

DIMENSION  B(4)»C(J. 

1  A*0*0 

2  DO  3  1  =  1  .4*1 

3  A=A+B ( I ) *C ( I ) 

RETURN 

END 

SUBROUTINE  MULTG ( A»B*C» I ) 

L  1ST 
LABEL 

this  subroutine  performs  the  dot  product  of  a  vector 
where  the  I  th  time  is  included 

DIMENSION  A  (  3  0  )  »B<30*4)*CU) 

1  A (  I  ) *0«  0 

2  DO  3  J*1  f4.1 

3  A  (  I  )*A< I )+B( I t J) *C( J ) 

RETURN 

END 

SUBROUTINE  MULTH(A*B.C) 

L  1ST 
LABEL 

this  subroutine  multiplies  a  matrix  by  a  row  vector 

To  GIVE  A  ROW  VECTOR 

1  D  0  4  I  *  1  *  4 , 1 

2  A  (  I )  =*0«0 

3  DO  4  J=1  *4*1 

4  A  (  I  )  =  A  (  I  )+B(J)*C  (  J » ]  ) 

RETURN 

END 

SUBROUTINE  MULTI  (AiBtCtl  ) 

L  1ST 
LABEL 

this  subroutine  multiplies  A  matrix  BY  A  ROW  VECTOR 
TO  GIVE  A  ROW  VECTOR  WHERE  THE  I  TH  TIME  IS  NEEDED 
DIMENSION  A(30»4)*B(30»4)»C(30»4t4) 

1  DO  4  J a  1  ♦  4 » 1 

2  A  (  I » vM  a  0  •  0 

3  DO  4  K=  1  *  4  » 1 

4  A ( ItJ)=A< I »J)+B< I »K)*C<I *K»J) 

return 

END 


n  n  *  *  n  *  *  on 


SUBROUTINE  SUMA < VX ♦ VY * VZ ) 


*  LIST 

LABEL 

THIS  SUBROUTINE  SUMS  ALL  OF  THE  MATRICES  INVOLVING 

the  i  th  time  and  divides  the  result  by  vz 

DIMENSION  VX  (  4  *  4  )  »  VY<30»4»4) 

5111  DO  5117  Js  1  *  4  *  1 

5112  DO  5117  K* 1  *  4  *  1 

5113  VX ( J  *K ) *0  *0 

5114  DO  5116  1*1*30*1 

5116  VX(J*K)  »  VX< J*K)+VY( I*J.K) 

5117  VX(J*K)*VX(J*K)/VZ 
RETURN 
END 


SUBROUTINE  SUMB(ZX.ZY.ZZ) 

L  1ST 
LABEL 

this  subroutine  sums  two  matrices 

DIMENSION  ZX(4*4)*ZY(4*4)*ZZ(4*4) 

7111  DO  7113  J* 1  *  4  *  1 

7112  DO  7113  K* 1  * 4  *  1 

7113  ZX( J»K)*ZY(j*K>+ZZ( J.K) 

return 

END 


SUBROUTINE  SUMC  < SX * SY * SZ * SW *  I ) 

LIST 

LA3EL 

this  subroutine  sums  two  matrices  involving  the 

I  TH  TIME  AND  DIVIDES  EACH  ELEMENT  BY  SW 
DIMENSION  SX( 30*4*4)  *SY( 30 *4*4) *SZ( 30 *4 *4) 

911  DO  913  J  =  1 *4  *  1 

912  DO  913  K=1 .4* 1 

913  SX< I  * J*K)  =  <SY<  I  * J*K)+SZ( I *J*K) )/SW 
RETURN 

END 


VJ  KJ 
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SUBROUTINE  TRACE<TA,TY) 

*  LIST 

*  LABEL 

C  THIS  SUBROUTINE  TAKES  THE  TRACE  OF  A  MATRIX 

DIMENSION  T  Y  (  4  #  4  ) 

6111  TA*0#0 

6112  DO  6113  J«1.4,l 

6113  TA*TA+TY( J»J) 

6114  PRINT  6 1  1 5  » T A 

6115  FORMAT  (  1H0*5X9H  TRACE  IS*E20.10) 

RETURN 

FND 


subroutine  traces  (Ta*ty»ij 

LIST 

LABEL 

this  subroutine  takes  the  trace  of  a  matrix  involving  the 

I  TH  TIME 

DIMENSION  T A ( 30) »TY( 30*4*4) 

781  TA*0« 

782  DO  786  J*1.4* 1 

783  Ta<  I  )cTA(  I  )+TY( I . J. J ) 

786  CONTINUE 

784  PRINT  78  5  »  TA( I) *  I 

785  FORMAT ( 1H0,5X9H  TRACE  I S * E20 . 1 0 t 5X 1 2H  TIME  NUMBER  *14) 
RETURN 

END 


non 
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SUBROUTINE  TRANSM  (  TM  ♦  F  I  *  I ) 

*  LIST 

*  LABEL 

c  this  subroutine  sets  up  the  state  transition  matrix 

DIMENSION  TM(30,4*4),T(30)*EI(4*4) 

COMMON  FA.W* SIGMAS, T ,EI • DELTMR . DELTML 
TM< I *1 *1 )«1.0+2.0»<SINF< (FA-FI )/2.0))**2 
TM (  1*1 *2)«SINF( FA-FI ) 

TM (  1*1 *3)«SINF( FA-FI )/W 

TM(  1*1  *4)* (4.0*(SINF( (FA-FI )/2*0)  )  **  2  )  /W 
T  M  (  I  *2  *  1  ) *— 3 • 0* (  FA-FI  )+2*0*SINF  (FA-FI  ) 

TMC 1*2 •2)*1*0-4#0»(SINF(  (FA-FI  )/2*0) )**2 
TM( I #2  *3 ) *-Tm( I *1*4) 

TM( I *2 *4) =-3,0* (FA-FI )/W+4*0*(SINF(FA-FI ) ) /W 
TM( 1*3*1 ) * 3 • 0*W* ( FA-F I ) -W*SlNF ( FA-FI ) 

TM( I*3*2)«2«0#W*(SlNF( (FA-FI )/2.0))**2.0 
TM( I *3  * 3 ) =TM ( 1*1*1) 

TM( 1*3*4) «-TM( 1*2*1) 

TM( I *4  *  1 ) ■  — TM ( 1*3*2) 

TM ( 1.4*2) *-W*TM ( 1*1.2) 

TM ( I .4*3) =-TM( I . 1.2 ) 

TM  (  I*4*4)*=TM(  1.2.2) 

RETURN 

END 


subroutine  TRANMI ( TMI *TM. I ) 

*  LIST 

*  LABEL 

THIS  SUBROUTINE  SETS  UP  THE  INVERSE  OF  THE  STATE  TRANSITION 
MATRIX  OR  THE  DERIVATIVE  OF  THE  INVERSE  OF  THE  STATE 
TRANSITION  MATRIX 

DIMENSION  TM( 30.4*4) *TMI ( 30*4.4) 

TMMI.1*1)*TM(I.3*3) 

TMI ( I , 1 .2) «TM( I .4*3 ) 

TMI ( I. 1*3) =  — TM (  1.1.3) 

TMI ( 1*1*4) x  — TM(  I  .2*3) 

TMI ( I  *  2  *  1 ) *TM( 1.3*4) 

TMI ( I  *  2*2) *TM( 1*4.4) 

TMI (1*2*3) «-TM(  1*1*4) 

TMI (1*2*4) *  — Tm (  1*2*4) 

TMI (1*3*1) «— TM ( 1*3*1) 

TMI ( I  *  3.2) ®-TM(  1.4,1) 

TMI ( I  *  3,3)«TM( I .1*1) 

TMI (  I  ,  3  *4 ) *  TM ( I ,2*1  ) 

TMI (1*4,1) *— TM ( 1*3,2) 

TMI ( I .4.2) *-TM( 1*4,2 ) 

TMI ( ! *4,3) =TM( 1*1,2) 

TMI ( I  *  4  *4 ) =TM ( 1*2*2) 

RETURN 
END 


u  *  *  o  u 
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SUBROUTINE  TRANIT  <TMIT*TM*I) 

*  LIST 

*  LABEL 

THIS  SUBROUTINE  SETS  Up  THE  TRANPOSED  INVERSE  OF  THE  STAT 

transition  matrix  or  the  derivative  of  the  tranpcsed  inve 
of  the  state  transition  matrix 

DIMENSION  TM  ( 30  *  4  *4 ) *  TM I T  (30*4*4) 


TM  I T  (  I 

♦  1 

.1 

)  *TM(  I  * 

3.3) 

TMlT  (  I 

♦  1 

*2 

)«TM(  I  * 

3.4  ) 

TM  I  T  (  I 

*1 

*3 

)»-TM(  I 

•  3.1  ) 

TMIT(  I 

*1 

*4 

)  *-TV(  I 

♦  3.2) 

TM  I T  (  I 

.2 

♦  1 

)  ■  T  M  (  I  * 

4.3) 

TM  IT  (  I 

♦  2 

.  0 

W  fc. 

)  «  TM  (  I  . 

4*4  ) 

TMlT  (  I 

♦  2 

♦  3 

)  »- TM (  1*4*1) 

T  M I  T  (  I 

*2 

*4 

)s-TM(  I 

*4.2) 

TMlT  (  I 

*3 

*1 

)  *  -  T  M  (  I 

*1*3) 

T  M  I T  (  I 

♦  3 

*2 

) «  —  TM (  I 

*1*4) 

TMlT  (  I 

♦  3 

*3 

)*TM(  I  * 

1*1) 

T  V  J  T  (  I 

*3 

*4 

)  =  tm(  i  * 

1*2) 

TMlT  (  I 

*4 

*2 

)  ■-  TM  (  I 

*2.3) 

TM !  T  (  I 

*4 

*2 

)  *-TM(  I 

*2.4) 

TM  I  T  (  I 

*4 

♦  3 

)  *  TM  (  I . 

2.1) 

TMI T  (  I 

*4 

*4 

)«TM(  I  * 

2*2) 

return 

END 

subroutine  dtranm  (DTM*Fl*n 

LIST 

LABEL 

this  subroutine  sets  up  the  derivative  of  each  element 
of  the  state  transition  matrix 

DIMENSION  DTM( 30.4*4) fEI (4*4) *T(30) 

COMMON  FA.W.SIGMAS.T.EI *  DEL TMR *  DEL TML 
DTM ( I  *  1  *  1 ) *-S I NF ( FA-F I ) 

DTM( I* 1  *  2 ) *-COSF( FA-FI ) 

DTM( I* l*3)«-(COSF(FA-FI >  ) /W 
DTM( I . 1 *4) *-2.*(S INF (FA-F I ) ) /W 
DTM( I  *  2* 1 ) *3*  — 2 • COSF (FA-FI ) 

DTM( I ♦ 2*2 ) =2»*SINF ( FA-FI ) 

DTM( I  *  2*3) r-DTM(  1.1*4) 

DTM ( I *2*4) =3*/W-4»*(COSF (FA-FI )  )  /W 
DTM(  I  *  3*  1 )  »-3.*W+W*COSF( FA-FI ) 

DTM ( I  *  3*2) «W*DTM(  1*1*1) 

DTM( I . 3.3)«DTM( 1*1.1) 

DTM( I ♦ 3*4) a-DTM ( 1.2.1) 

DTM( I *4*1)«W*SINF (FA-FI ) 

DTM( I  *  4*  2 ) =W*COSF ( FA-FI ) 

DTM ( I*4*3)cCOSF(FA-FI ) 

D  TM ( I *4*4)«DTM( 1*2*2) 

RETURN 

END 


SUBROUTINE  AOM(OM) 


*  LIST 

*  LABEL 

C  THIS  SUBROUTINE  SETS  UP  THE  Q  MATRIX 

DIMENSION  QM ( 4  »  4  ) 

DO  3  1=1 *4*1 
DO  3  Js 1 *4*1 
3  OM ( I  *  J ) *0.0 
OM (1*1) *1.0 
QM(2*2  )  «*1.0 
RETURN 
END 


LCD 


SUBROUTINE  ERROR  I ( El  E ) 

*  L  1ST 

*  LABEL 

C  THIS  SUBROUTINE  TRANSFERS  THE  INITIAL  ERROR  TO  THE 

C  TARGET  CO-ORDINATE  FRAME 

DIMENSION  E I  E  (  4  *  4  ) *  TM ( 4  *  4 )  *  TM  T ( 4  *  4 ) • A  <  4  *  4 ) » 

1EI (4*4 ) *  T ( 3 0  ) 

COMMON  FA*W*SIGMAS.T  tEI *  DEL TMR * DELTML 
TM(ltl  )**  1.0  +  2.0  <  S  INF  <  FA/2.0)  )  **2 
TM(1*2)=SINF(FA) 

TM( 1 *3 ) =SINF( FA)  /W 

TM(  1  *4 ) *4.0* ( SINF( FA/2.0 )**2 ) /W 

TM ( 2  *  1 ) *-3.0*fa  +  2.0*SINF(FA) 

TM ( 2  *  2 ) =1.0-4.0*( SINP( FA/2.0 )  ) **2 
TM  (  2  *  3  ) *-TM( 1 *4) 

TM ( 2  *  4 ) *(-3.0*FA  +  4.0*SlNF(  FA)  ) /W 
TM ( 3 1 1 )«3.0*W*FA-W*SINF(FA) 
TM(3*2)«2.0*W*(SlNF(FA/2.0) )**2 
TM(3*3)STM(1*1) 

TM ( 3  *  4 ) «— TM (2*1) 

T  M ( 4  *  1 ) *-TM( 3*2 ) 

TM(4*2 ) *-W*TM (1*2) 

TM ( 4  *  3 ) *-TM( 1*2 ) 

TM ( 4  *  4 ) *  TM (2.2) 

TMT(1,1)«TM(1.1) 

TMT(1,2)*TM(2.1) 

TMT ( 1 .3 ) *TM( 3* 1 ) 

TMT ( 1  *  4 ) *  TM ( 4 . 1 ) 

TMT(2*1)«TM(1»2) 

TMT ( 2  *  2 ) *TM( 2*2 ) 

TMT(2*3 )*TM( 3*2 ) 

TMT (2*4)= TM (4*2) 

TMT  (  3  *  1  )  *=TM  (  1*3) 

TMT(3*2)=TM(2.3) 

TMT(3.3)=TM(3.3) 

TMT ( 3*4 ) =TM( 4*  3 ) 

TMT ( 4 » 1 ) =TM( 1*4) 

TMT ( 4  *  2 ) =TM  (  2*4) 

TMT(4.3)*TM(3»4) 

TMT (4*4 ) =TM( 4*4) 

CALL  MULTB( A.EI .TMT) 

CALL  MULTB(EIE*TM*A) 

RETURN 

END 


n  n 


SUBROUTINE  ERROR M ( EA* H . COST  * T T ) 


C  BLOCK  ONE 


♦  LIST 


LABEL 

this  subroutine  computes  the  target  estimation 

ERROR  COVARIANCE  MATRIX  AND  THE  COST 

DIMENSION  C(30*4*4)»CM(4*4)*E(30*4*4)*EA(4*4) * EAC (4*4) 

1  EAI(4.4)»EI(4»4)»EII(4»4).EIE(4.4).H(4)*HH(4.4)* 

2  TM(30*4*4)*TMI(30*4*4)*TMIT(30*4*4)*T(  30)*TT(30) 
COMMON  F  A  *  W ♦ S I GMAS  »  T ♦ E I *DELTMR  * DELTML 

810  CALL  ERROR  I ( EIE) 

816  CALL  MULTA  (HH*H) 

811  DO  818  1=1 *30*1 

812  F I »W*TT {  I  ) 

813  CALL  TRANSM  (TM*FI.I) 

814  CALL  TRANMI  {TMI*TM*I) 

815  CALL  TRANIT  (TMIT.TM*I) 

817  CALL  MULTC  (C*HH*TMI*I) 

818  CALL  MULTD  (E*TMlT*C*I> 

819  CALL  SUMA  ( EAC*  E*S IGMAS ) 

8191  DO  8194  I ■ 1  *  4  *  1 

8192  DO  8194  J=l*4*l 

8193  El  I ( I  * J)  *  0.0 

8194  E I  I ( I  *  J ) =E I E (  I*J) 

820  CALL  INVERT  ( 4 ♦ E 1 1 ) 

821  CALL  SUMB  ( EA I  *  E I  I  *  EAC  > 

8211  DO  8214  I  *=  1  *4*1 

8212  DO  8214  J-l*4*l 

8213  EA ( I  *  J )  *  0.0 

8214  EA ( I  * J )  *  EA I ( I  .  J  ) 

822  CALL  INVERT  (4*EA) 

823  CALL  AOM(OM) 

824  CALL  MULTB(CM*OM.EA) 

825  CALL  TRACE ( COST  *  CM ) 

RETURN 

END 


r»  o 


SUBROUTINE  ERRORM  (EA»H*COST> 


C  BLOCK  TWO 


# 

# 


810 

811 
812 
813 
8 1 A 

815 

816 

817 

818 

819 

8191 

8192 

8193 

8194 

820 
821 

8211 

8212 

8213 

8214 
822 

823 

824 

825 


LIST 

LABEL 

this  subroutine  computes  the  target  estimation 

ERROR  COVARIANCE  MATRIX  AND  THE  COST 

DIMENSION  C(3C»4,4) ,CM(4,4) ,E(30,4,4) ,EA(4,4) ,EAC(4,4) 

1  EAI ( 4  *  4 ) *  E I (4*4)*EIE(4,4) »  E  I I (4*4) ,  H  (  3  0  (  4  ) , HH (30,4,4) 

2  TM(30,4,4),TMI (30,4,4) ,TMlT( 30,4,4) ,T(30) 

COMMON  FA,W,SIGMAS,T,EI .DELTMR.DELTML 
CALL  ERROR  I ( E I E ) 

DO  818  1  =  1  .30,1 
F I =W*T ( I ) 

CALL  TRANSM  (TM.FI, I) 

CALL  TRANMI  ( TMI ,TM* 1 ) 

CALL  TRANIT  (TMIT.TM, I ) 

CALL  MULTAAtHH.H, I ) 

CALL  MULTD(C,HH,TMI , I ) 

CALL  MULTO  ( E  »TM I T , C , I ) 

CALL  SUMA  (EAOE.SIGMAS) 

DO  8194  1*1.4, 1 
DO  8194  J  =  1 » 4  » 1 
EII(I.J)  =  0.0 
EII(I,J)=EIE(  I.J) 

CALL  INVERT  (4, Eli) 

CALL  SUMB  (EAI.EII ,EAC) 

DO  8214  1*1, 4.1 
DO  8214  J* 1 , 4 , 1 
E A ( I ♦ J )  =  0,0 
EA(I.J)  *  EAKI.J) 

CALL  INVERT  ( 4 *  E A ) 

CALL  AOM(OM) 

CALL  MULTB(CM.OM.EA) 

CALL  TRACE (COST.CM) 

RETURN 

END 


SUBROUTINE  HOR IZS( IHV.HR *HL *HS ) 

LIST 

LABEL 

DIMENSION  T(  3O ) *EI  (4*4 ) *  IHV( 30 ) ,HR(4) »HL<  4) *HS( 30*4  ) 

COMMON  FA*W*SIGMAS*T  *  El *DELTMR.DELTML 

THIS  SUBROUTINE  SETS  UP  THE  HORIZON  VECTOR  SCHEDULE 

1  DO  10  1=1*30*1 

2  IF  (IHV(I)-l)  3*7*7 

3  DO  5  Js 1  *  4  » 1 

4  HS ( I *J) *0.0 

5  H S ( I  *  J ) =HL ( J ) 

6  GO  TO  10 

7  DO  9  J=  1  *  4  *  1 

8  HS ( I*J)=0.0 

9  HS ( I  *  J ) =HR ( J ) 

10  CONTINUE 

RETURN 

END 


SUBROUTINE  LAMDA ( TLAMDA *HS *  I ) 

LIST 

LABEL 

this  subroutine  computes  the  lamda  vector  which  is  used  to 

COMPUTE  THE  INFLUENCE  COEF.  DUE  TO  A  CHANGE  IN  HORIZON 
DIMENSION  Tm( 30*4*4) *TMI (30*4*4) *TMlT(30*4*4) * TLAMDA ( 30  *  4 ) 
1EA(4»4)*QM(4*4)*T(30)»EI(4*4)*HS(30*4)*A(30*4*4)*B(30*4*4) 
2C ( 30  *  4 . 4 ) *D( 30*4*4) 

COMMON  FA *W* SIGMAS, T  *EI *  DEL TMR * DELTML 

4  CALL  ERRORM(EA*HS*COST) 

89  F I *W*T (  I  ) 

1  CALL  TRANSM ( TM  *  F  I  » I ) 

2  CALL  TRANM I ( TM  I  *  TM  *  I  ) 

3  CALL  TRANIT(TMIT,TM* I ) 

5  CALL  AQM(QM) 

6  CALL  MULTC( A*EA»TMIT*  I  ) 

7  CALL  MULTC(B*OM* A*  I  ) 

8  CALL  MULTC(C.EA*B* I ) 

9  CALL  MULTD(D*TMI*C*I  ) 

10  CALL  MULTI  ( TLAMDA*HS*D* I ) 

11  DO  12  J* 1  *  4  *  1 

12  TLAMDA (I  *  J  )*- 2 • 0* TLAMDA ( I  *  J ) /S IGMAS 
RETURN 

END 


c 


SUBROUTINE  DERROR  ( ERR  »H  *  I  *  TT  ) 
BLOCK  ONE 


* 

* 


981 

982 

983 

984 

985 

986 

987 

988 

989 

990 

991 

992 

993 


LIST 

LABFL 

DIMENSION  A(30*4*4)*B<30*4*4)*D(30*4*4) *DTM(30*4*4) * 

1  DTMI ( 30*4*4) »DTMIT ( 30*4*4) *  E  I  (  4  *4  ) *  ERR ( 30*44 ) » 

2  F(30.4*4)*T(30).TM(30*4*4)*TmK30*4*4)*TMIT<30*4*4) 

3  H ( 4 ) ,HH( 4*4) *TT ( 30 ) 

COMMON  FA  *  W*  S I GMAS  *  T *EI * DELT MR  *  DEL TML 
F I =W*T  T ( I ) 

CALL  DTRANM(DTM*FI *  I ) 

CALL  TRANMI  (DTMI*DTM*I) 

CALL  TRANIT  (DTMIT*DTM.I ) 

CALL  TRANSM  (  Tm  *  F I  *  T ) 

CALL  TRANMI ( TMI *  TM* I ) 

CALL  TRANIT  ( TM I T  *  TM  *  I ) 

CALL  MULTA  (HH*H) 

CALL  MULTC( A*HH*TMI , I ) 

CALL  MULTC  ( B*HH  *  DTM I  *  I ) 

CALL  MULTD  ( D*DTM IT , A* I ) 

CALL  MULTD  ( F*TMIT*B* I ) 

CALL  SUMC  (ERR.D*F*SIGMAS*I ) 

RETURN 

END 


130 


SUBROUTINE  DERROR ( ERR , H *  I) 

C  BLOCK  TWO 

*  LIST 

*  LABEL 

DIMENSION  A  ( 30*4*4) »  8  ( 30*4*4) *0(30*4*4) *DTM (30*4*4), 

1  DTM I ( 30*4*4) *DTMlT( 30,4*4) *EI (4*4 ) *  ERR ( 30*44 ) » 

2  F(30,4*4)*T{30)*TM{30*4.4)*TMI(30.4*4).TMIT(30*4*4), 

3  K( 30,4 ) ,HH( 30*4*4) 

COMMON  FA*W,SIGMAS,T ,EI » DELTMR *DELTML 

981  F I =W#T ( I ) 

982  CALL  DTRANM(DTM*FI , I  ) 

983  CALL  TRANMI  (DTMI*DTM#I) 

984  CALL  TRANIT  ( DTM I T * DTM *  I  ) 

985  CALL  TRANFM  (TM,FI*I) 

986  CALL  TRANMI ( TMI »TM, I  ) 

987  CALL  TRANIT  (TMIT*TM*I) 

988  CALL  MULT AA ( HH*H* I ) 

989  CALL  MUL  TD ( A*  HH  *  TM I  *  I  ) 

990  CALL  MULTD( B*HH* DTMI *  I ) 

991  CALL  MULTD  ( D *DTMI T , A  *  I  ) 

992  CALL  MULTD  (F*TMIT*B*I) 

993  CALL  SUMC  ( ERR  * D * F * S I GMAS , I ) 

RETURN 

END 


r»  n  r> 
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SUBROUTINE  CHECK  ( T AR I M *  TAR  I * OLCOSR , HR * T T , CRUD ) 

C  BLOCK  ONE 

CHANGING  ONLY  ONE  TIME  SHOULD  RESULT  IA  AN 
ACCURATE  PREDICTION  OF  CHANGE  IN  COST 
WITHIN  .10/0  OF  ACTUAL  CHANGE 
L  1ST 
LABEL 

DIMENSION  EAR (4. 4) *EI(4,4)*F(30), 

1  TAR  I  (  30 ) *  T ( 30 ) 

COMMON  FA*W,SIGMAS,T,EI *DELTMR  *DELTML 
PRINT  16 

16  FORMAT ( 1H5 . 10X14H  ACCURACY  TEST) 

1  B IGA=ABSF( TARIM) /CRUD 

2  TT(15)=TT( 15J-TARI (155/BIGA 

3  CALL  ERRORM  ( EAR *HR * PCOSR , TT ) 

CALL  PRINTS  (EAR) 

4  COSTR=(TARl ( 15)*** )/BIGA 
PRINT  14*  COS  TR 

14  FORMAT ( 1H0 * 10X28H  PREDICTED  COST  CHANGE  TEST,E20.8) 

5  ACOSTRaOLCOSR-P :osr 
PRINT  15*  ACOSTR 

15  FORMAT( 1H0, 10X25H  ACTUAL  COST  CHANGE  TEST*E20.8) 

6  I F ( ACOSTR )  7*7*8 

7  PRINT  71 

71  FORMAT( 1H0.10X38H  NEW  COST  GREATER  THAN  OLD  COST  TEST) 
GO  TO  131 

0  I F ( (COSTR- ACOSTR ) /OLCOSR )  81*13*91 
81  BUTT= ( ACOSTr-COSTR ) /OLCOSR 
9  IF  (BUTT-.001)  13*13,11 
91  RUTT= ( COSTR-ACOSTR ) /OLCOSR 

10  IF  ( RUTT- #  00 1 )  13*13,12 

11  PRINT  111, BUTT 

111  FORMAT ( 1H0 , 10X32H  TEST  ACTUAL  EXCEEDS  PREDICTED* 

1  KX ,F20. 10 , 5X19H  PERCENT  DIFFERENCE) 

GO  TO  131 

12  PRINT  1 2  1  • RUT  T 

121  FORMAT ( 1H0, 10X32H  TEST  PREDICTED  EXCEEDS  ACTUAL* 

1  5X*F20.10,5X19H  PERCENT  DIFFERENCE) 

13  PRINT  132 

132  FORMAT ( 1H0, 10X27H  SO  FAR  SO  GOOD  TEST  WORKS) 

131  TT( 15) =TT( 15) +TARI ( 15 J/BIGA 
RETURN 
END 


132 


SUBROUTINE  CHECK  < T AR I M *  TAR  I .OLCOSR . HS ) 

C  BLOCK  TWO 

C  CHANGING  ONLY  ONE  TIME  SHOULD  RESULT  IA  AN 

C  ACCURATE  PREDICTION  OF  CHANGE  IN  COST 

C  WITHIN  .10/0  OF  ACTUAL  CHANGE 

*  LIST 

*  LABEL 

DIMENSION  TARM  30  )  t  F<  30)  •  T(  30)  tEI  Ut4)  *EA<4t4  )  tHS(  30#4) 
COMMON  FA.W.SIGMAS.T.EI .DELTMR . DELTML 
PRINT  16 

16  FORMAT( 1H5 .10X14H  ACCURACY  TEST) 

1  8IGA*ABSF( TARIM) /DELTMR 

2  T ( 1 5 ) *T { 15) “TAR  I ( 1 5 ) /B I GA 

3  CALL  ERRORM  ( EA * HS * PCOSR  ) 

CALL  PRINTS  ( EA ) 

4  COSTr=(TARI ( 15)**2)/BIGA 
PRINT  1 4 . COSTR 

14  F0RMAT(1H0*10X28H  PREDICTED  COST  CHANGE  TEST.E20.8) 

5  ACOSTr*OLCOSR-PCOSR 
PRINT  15*  ACOSTR 

15  FORMAT ( 1H0* 10X25H  ACTUAL  COST  CHANGE  TEST.E20.8) 

6  I F ( ACOSTR )  7.7.8 

7  PRINT  71 

71  FORMAT (1H0*10X38H  NEW  COST  GREATER  THAN  OLD  COST  TEST) 
GO  TO  131 

8  I F ( (COSTR- ACOSTR) /OLCOSR )  81.13.91 
81  BUTT=( ACOSTR- COSTR ) /OLCOSR 

9  IF  (BUTT-. 001)  13.13.11 

91  RUTT= ( COSTR-ACOSTR ) /OLCOSR 

10  IF  (RUTT-.OOl)  13.13.12 

11  PRINT  111. BUTT 

111  FORMAT( 1H0. 1CX32H  TEST  ACTUAL  EXCEEDS  PREDICTED. 

2  5X.F20.10.5X19H  PERCENT  DIFFERENCE) 

GO  TO  131 

12  PRINT  1 2 1 . RUT  T 

121  FORMAT(1HO,10X32H  TEST  PREDICTED  EXCEEDS  ACTUAL. 

1  5X.F20.10.5X19H  PERCENT  DIFFERENCE) 

13  PRINT  132 

132  FORMAT( 1H0, 10X27H  SO  FAR  SO  GOOD  TEST  WORKS) 

131  T( 15)*T( 15 J+TARI < 15 )/BlGA 
RETURN 
END 


*  n  r> 
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SUBROUTINE  LOG  I C ( B I GA *  TAR  I M * CRUD* T AR I *OLCOSR , PCOSR * 
1  EAR  *HR  *  TT  ) 

C  BLOCK  ONE 

determines  new  measurement  time  schedule  AND 

NEW  COST 
L  1ST 
LABFL 

DIMENSION  F  (  3  0  )  t  EAR  (4*4)  *  E  I  <  4  *  4  )  *  HR  (  4  )  *  FR  (  30  )  » 

1  DELTR(30) *T(30)*TR(30)*TT(30) 

COMMON  FA.W*SIGMAS*T*EI * DELTMR * DELTML 

1  BIGA  =  ABSF(TARIM)  /CRUD 
PRINT  26  *  CRUD 

26  FORMAT ( 1H0,10X19H  MAX  TIME  I N CREMENT , 5X , F20 . 1 0 ) 
PRINT  27  »  B I GA 

27  FORMAT ( 1HO,10X13H  SCALE  FACTOR  * 5X * E20*8 ) 

2  DO  45  1*1*30,1 

3  DELTp (  I )  =  T AR I ( I )/BIGA 

4  T R ( I )  =  T T ( I)-DELTR(I) 

41  I F ( TR (  I )) 42  *  45  *  43 

42  DELTr ( I ) *DELTR ( I ) +Tr ( I ) 

TR(I)«0.0 
GO  TO  45 

43  I F ( TR ( I ) - ( FA/W ) )  45*45*44 

44  DELTR  (  I  UDELTR(  I  )+ ( TR ( I ) -( FA/W  )  ) 

TR ( I )=FA/W 

45  CONTINUE 

5  DO  6  1  =  1  *30*1 

61  F ( I ) =W*  TR ( I ) 

62  F(I)*F( I)*180./3. 14159 

6  T  T ( I ) *TR ( I ) 

7  CALL  ERRORM  ( EAR *HR , PCOSR * TT ) 

8  C0STR=0 • 

9  DO  10  1=1*30,1 

10  COSTRsCOSTR  +  TARI ( I )»DELTR(  I  ) 

11  ACOSTR -OL COS R -PCOSR 
PRINT  2  8, COS T R 

28  FORMAT ( 1HO,10X22H  PREDICTED  COST  CHANGE  * 5X , E20. 8 ) 
PRINT  29  *  ACOSTR 

29  FORMATt 1H0,10X19H  ACTUAL  COST  CHANGE  * 5X , E20 . 8 ) 

12  IF  (COSTR-ACOSTR)  131*15*141 
131  BOO=ACOSTR/COSTR 

13  IF  ( BOO- 10. )  15*15*19 
141  FOO=COSTR/ACOSTR 

14  IF  (FOO-lO.)  15*15*20 

15  APERR=ACOSTR/OLCOSR 
PRINT  30  *  APERR 

30  FORMAT ( 1H0*10X15H  PERCENT  CHANGE  * 5X , F20 . 10 ) 

IF  (APERR) 111 *333*555 

111  IF  ( .001+APERR ) 222*333*444 


222  CRUD=CRUD*< *00l/ABSF( APERR) ) 

GO  TO  1 

333  CRUD=CRUD/2*0 
GO  TO  1 

444  CRUD*CRUD*( 1.0-ABSF( APERR) /•OOl ) 

GO  TO  1 

555  IF  (APERR-. 001)21*18*18 

21  CRUD=CRUD* ( •001/APERR ) 

IF(CRUD-10000. ) 22 *18 *18 

22  GO  TO  1 

18  PRINT  31 

31  FORMAT! 1H0.10X18H  NEW  TIME  SCHEDULE) 

CALL  PRINTA(TT) 

PUNCH  3 1 1  *  TT 

311  FORMAT ( F20  # 1 0 ) 

PRINT  32 

32  FORMAT  (  1H0»10X19H  NEW  CENTRAL  ANGLES) 

181  CALL  PRINTA  (F) 

PRINT  33 

33  FORMAT! 1H0,10X19H  FINAL  ERROR  MATRIX) 

CALL  PRINTS  (EAR) 

182  CALL  PRINTR  ( OLCOSR ,PCOSR *COSTR * ACOSTR ) 

GO  TO  25 

19  PRINT  191*  BOO 

191  FORMAT ( 1  HO  * 10X40H  ACTUAL  GREATER  THAN  TEN  TIMES  PREDICTED 
1  F20.10) 

GO  TO  25 

20  PRINT  201  *  FOO 

201  FORMAT!  1H0*10X40H  PREDICTED  GREATER  THAN  TEN  TIMES  ACTUAL 
1  F20.10) 

25  CONTINUE 
RETURN 
END 


♦  on 
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ED* 


AL* 


SUBROUT INE  logici biga *tarim*crud*tari *OLCOSR *pcosr*ea*h1 
C  BLOCK  TWO  I 


DETERMINES  NEW  MEASUREMENT  TIME  SCHEDULE  AND 
NEW  COST 
L  1ST 
*  LABEL 

DIMENSION  F(30)*TR(30)»FR(30) *DELTR ( 30)  .TARI ( 30  )  * 
1EA(4*4)*HS(30*4)*T(30)*EI(4*4) 

COMMON  FA*W*SIGMAS»T  *EI *DELTMR *DELTML 

1  B  IGA*A8SF ( TARIM)/ CRUD 
PRINT  2  6  *  CRUD 

26  FORMAT! 1H0, 10X19H  MAX  TIME  I NCREMENT * 5X * F20 . 1 0 ) 
PRINT  27  * B I GA 

27  FORMAT! 1H0*10X13H  SCALE  FACTOR  * 5X , E20 . 8 ) 

2  DO  45  1*1*30,1 

3  DELTR! I ) =TARI ( I ) /BIGA 

4  T R ( I ) *T  (  I ) -DELTR  < I ) 

41  I F ( TR {  I  )  ) 42  *  45  »  43 

42  DELTR (  I ) *DELTR( I )+Tr !  I  ) 

TR! I )=0*0 

GO  TO  45 

43  I F ( TR ( I ) -( FA/W) )  45*45*44 

44  DELTR ! I ) *DFLTr ( I ) +( T R ( I ) - ( FA/W ) ) 

TR! I )=FA/W 

45  CONTINUE 

5  DO  6  1=1 *30,1 

61  F ( I ) =W* TR (  I  ) 

62  F(I ) * F ( I  )*180#/3. 14159 

6  T ( I  )  *TR (  I  ) 

7  CALL  ERRORM  ( EA *HS , PCOSR  ) 

8  COSTR  =  C  • 

9  DO  10  1=1*30*1 

10  COSTR=COSTR  +  TARI ! I )* DELTR!  I  ) 

11  acostr^olcosr-pcosr 

PRINT  2 8  * COSTR 

28  FORMAT ( 1H0*10X22H  PREDICTED  COST  CHANGE  * 5X , E2 0* 8 ) 
PRINT  29 , ACOSTR 

29  FORMAT! 1H0*10X19H  ACTUAL  COST  CHANGE  * 5X , E 20 . 8  ) 

12  IF  (COSTR-ACOSTR)  131*15,141 
131  BOO=ACOSTR/COSTR 

13  IF  (BOO-IO.)  15*15*19 
141  F  00=C0S  T  R/ACOSTR 

14  IF  (FOO-lO.)  15*15*20 

15  APERR=ACOSTR/OLCOSR 
PRINT  30  ,  APERR 

30  PORMAT! 1H0*10X15H  PERCENT  CHANGE  * 5X . F20 ♦ 1 0 ) 

I  F ( APERR )111*333*555 

111  IF!  •  001  -t-AP  ERR)  2  22  *333  *444 
222  CRUD=CRUD*( .OOl/ABSF!  APERR)  ) 


GO  TO  1 

333  CRUD=CRUD/2 • 0 
GO  TO  1 

444  CRUD=CRUD*Q«0-ABSF(  APSRR  >  /  •  00 1  ) 

GO  TO  1 

555  IF( APERR~#OOl )21 ,18,18 

21  CRUD=CRUD*(.OOl/APERR) 

IF(CRUD-10000.  )22#18tl8 

22  GO  TO  1 

18  PRINT  31 

31  FORMAT ( 1HO*10X18H  NEW  TIME  SCHEDULE) 

CALL  PRINTA  (T) 

PUNCH  3 1 1 »  T 
311  FORMAT ( F20# 10 ) 

PRINT  32 

32  FORMAT ( 1H0* 10X19H  NEW  CENTRAL  ANGLES) 

181  CALL  PRINTA  (F) 

PRINT  33 

33  FORMAT ( 1HO»10X19H  FINAL  ERROR  MATRIX) 

CALL  PRINTS  ( EA ) 

182  CALL  PRINTR  (  OLCOSR , PCOSR * COSTR » ACOSTR ) 

GO  TO  25 

19  PRINT  191* BOO  - 

191  FORMAT ( 1H0»iOX4OH  ACTUAL  GREATER  THAN  TEN  TIMES  PREDICTED 
1  F20.10) 

GO  TO  25 

20  PRINT  201, FOO 

20 1  FORMAT ( 1HO»10X40H  PREDICTED  GREATER  THAN  TEN  TIMES  ACTUAL 
1  F20.10) 

25  CONTINUE 
RETURN 
END 
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SUBROUTINE  jump  ( TLAMDA*DHLR*DHRL»PCOSR , IHV.HR.HL ) 

*  LIST 

*  LABFL 

DIMENSION  IHV( 30 ) .DELCOS ( 30) * TLAMDA ( 30 , 4 ) . EA ( 4 ,4 ) , 

1HS ( 30,4  ) ,  DHRL ( 4 ) * DHLR ( 4 ) , HR(4),HL(4),T ( 30 ) , E I ( 4 ,4  ) 

COMMON  FA.W.SIGMAS.T  , El .DELTMR.DELTML 

PRED=0»  0 

SCALF=0.0 

NUM=0 

CALL  HOR IZS( IHV.HR.HL »HS ) 

DO  6  1=1*30*1 

CALL  LAMDA (TLAMDA »HSv I ) 

IFUHV(I)-l)  1.2*2 

1  CALL  MULTG  (DELCOS. TLAMDA *DHLR *  I ) 

GO  TO  3 

2  CALL  MULTG  ( DELCOS .TLAMDA .DHRL . I ) 

3  PRINT  81  *  I  *  DELCOSt  I ) 

81  FORMAT ( 1HO,10X20H  COST  CHANGE  FOR  THE*I4,6H  POINT .E20. 10 ) 
I F ( DELCOS (  I) ) 4 . 6  *  6 

4  IF(DELCOS(  I)-SCALE)5»6»6 

5  SCALE  =  DEL COS ( I ) 

NUM  =  I 

6  CONTINUE 

I  F  (  IHV  (  NL'M  )-l  >7.8*8 

7  I HV ( NUM ) = 1 
GO  TO  9 

8  I HV ( NUM ) =0 

9  PRED*-DELCOS(NUM) 

PRINT  100.PRED 

100  FORMAT ( 1H0.10X22H  PREDICTED  COST  CHANGE . E20 . 1 0 ) 

CALL  PR  I  NT  D  (  IHV ) 

CALL  HORIZS  ( IHV. HR , HL.HS ) 

CALL  ERRORM  (EA.HS.HOR) 

ACT  =  PCOSR— HOR 
PRINT  91 

91  FORMAT  (1H5.10X17H  NEW  ERROR  MATRIX) 

CALL  PRINTB(EA) 

IF  (ACT)  10.10,12 

10  PRINT  11 .ACT 

11  FORMAT ( 1HO.10X20H  SOMETHING  IS  WRONG  . 

1  24H  THE  NEW  COST  IS  GREATER, 5X.E20. 10) 

1E20.10) 

GO  TO  19 

12  PRINT  200, ACT 

200  FORMAT! 1H0.10X19H  ACTUAL  COST  CHANGE » E20. 10 ) 

DEL  =PRED/ ACT 
PRINT  14. DEL 

14  FORMAT ( 1HO  .10X29H  RATIO  OF  PREDICTED  TO  ACTUAL » F20. 10 ) 

19  CONTINUE 
RETURN 
END 


SUBROUTINE  printa(aa) 

*  LIST 

*  LABEL 

C  THIS  SUBROUTINE  PRINTS  OUT  A  LIST 

DIMENSION  AA { 3 0 ) 

PRINT  2  *  AA 
2  FORMAT ( 4E30 • 1 0 ) 

RETURN 

END 


SUBROUTINE  PRINTB(BB) 

*  LIST 

*  LABEL 

c  this  subroutine  prints  out  a  matrix 

DIMENSION  BB (4*4) 

PRINT  3 

3  FORMAT! 1  HO  * 10H  I  J  .20X10H  I  J  . 

1  20X10H  I  J  ♦  20X 1  OH  I  J  ) 

L  3 1 
LL  =  2 
LLL  =  3 
LLLL=4 

DO  5  I  =  1  *  4  *  1 

PRINT  4* I »LfBB(  I  »L ) f I vLL#BB( I *LL) f I  *  LLL  *  BB ( I »LLL ) ♦ 
ll.LLLLfBB! I *LLLL ) 

4  FORMAT (1HO*4(I5*I5*E20#10)  ) 

5  CONTINUE 
RETURN 
END 


SUBROUTINE  PRINTC(BBtl) 

*  LIST 

*  LABEL 


this  subroutine 

PRINTS 

OUT 

ONE 

OF  THE 

I  TH  MATR  ICES 

DIMENSION  BB ( 30 
PRINT  4 

*4  f  4  ) 

4  FORMAT! 1H0 * 5H 

I  *  5H 

J 

» 5  H 

K  * 

1  15X5H  I  *  5H 

J  1 5H 

K 

♦  1  5X5H  I 

*  5H  J  *  5H  K 

2  15X5H  I  *  5H 

J  ♦  5H 

K 

) 

L  3 1 

LL  =  2 

LLLa3 

LLLL34 

DO  8  J31 *4»1 

PPINT  7t I t JtL*BBt I  * J*L  >  1 1 1 J*LL*BB( I • JtLL) » I *J*LLLt 
2  BB ( I  *  J  »  LLL ) f I ♦ J  >LLLL  *  BB ( I ♦ J ♦ LLLL ) 

7  F0RMAT(1H0»4(3I5»E15.8)  ) 

8  CONTINUE 
RETURN 
END 


n  o  o 
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SUBROUTINE  PRINTDUHV) 

LIST 

LABEL 

iHV(n  is  the  horizon  vector  schedule 

IF  IHV(I)  *1  WE  ARE  USING  THE  RIGHT  HORIZON 
IF  IHV(I)  =0  WE  ARE  USING  THE  LEFT  HORIZON 
DIMENSION  I HV ( 30 ) 

DO  717  1=1,30,1 

711  IFUHV(I)-l)  712*715*715 

712  PRINT  713*1 

T13  FORMAT ( 1H0 , 5X7H  AT  THE*I3*25H  MEASUREMENT  POINT  WE  ARE* 
1  22 H  USING  THE  LEFT  HORIZON) 

714  GO  TO  717 

715  PRINT  716.1 

716  FORMAT ( 1H0.5X7H  AT  ThE*I3*25H  MEASUREMENT  POINT  WE  ARE* 
1  23H  USING  THE  RIGHT  HORIZON) 

717  CONTINUE 
RETURN 
END 


SUBROUTINE  PRINTR(A.B*C*D) 

*  LIST 

*  LABEL 
PRINT  1 

1  FORMAT ( 1H0,10X9H  OLD  COST,21X9H  NEW  COST, 
1  2 1X12H  PRED  CHANGE, 18X14H  ACTUAL  CHANGE) 

PRINT  2*A,B,C*D 

2  FORMAT ( 4E30 # 1 0 ) 

RETURN 

END 
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